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MSC/NASTRAN's performance in the solution of fully-coupled fluid/structure problems is 
evaluated. NASTRAN is used to perform normal modes (SOL 103) and forced-response analyses 
(SOL 108, 111) on cylindrical and cubic fluid/structure models. Each model is discretized using 
finite element methods. Bulk data file cards unique to the specification of a fluid/structure model 
are discussed and analytic partially-coupled solutions are derived for each type of problem. 
These solutions are used to evaluate NASTRAN's solutions for accuracy. Appendices to this 
work include NASTRAN data presented in fringe plot form, FORTRAN source code listings 
written in support of this work, and NASTRAN data file usage requirements for each analysis. 
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Length of a side for cube, plate 
Force amplitude 
Young’s Modulus 
Acoustic pressure 
Radius for cylinder 

Acoustic speed of sound in an ideal gas 
Thickness of plate, shell 

v^T 

Length of cylinder 
Time 

In-plane displacements 
Out-of-plane displacements 
Rectangular coordinate system 
Cylindrical coordinate system 
Damping coefficient 
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Section 1: Introduction 

With the advent of more powerful computers, it is becoming possible to solve vibrations 
problems of ever increasing complexity. A greater dependence upon computers necessitates 
that researchers know the limits and capabilities of the computer systems and the software 
packages they are using. If these limits are unknown, the researcher may be unable to discern 
nonsensical data from accurate data provided by the computer. 

It is the aim of this paper to evaluate the capabilities of a commercially available finite 
element analysis package called MSC/NASTRAN (versions 67 and 68) in the analysis of fully 
coupled fluid/structure problems. Normal modes analysis and forced response analysis is used 
on two simple, well-understood geometric models. The models are different from each other 
both in shape and complexity in an effort to determine the accuracy of NASTRAN. The 
accuracy of a NASTRAN-generated solution is determined through comparison with classical 
analytic solutions. 

This paper is divided into two sections, one each for a cubic and a cylindrical geometry. 
Within each of these sections are analytic and numeric solutions for three types of problems. 
The first problem is a normal modes analysis for a model containing fluid elements only. The 
second is a normal modes analysis of a model containing both fluid and structure elements. 
Normal modes of vibration are determined for the structural and fluid portions of the system. 
The third problem is a forced response analysis of the same model as that used in the second 
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problem. For the normal modes problems, both quadratic and linear element meshes are used. 
The forced response analysis makes use of a linear finite element mesh only. 

MSC/PATRAN was used as a graphical pre- and post-processor. FEM models and bulk 
data files were constructed within PATRAN. However, at the time of this work, PATRAN did 
not have the capability to work with fluid elements. Thus, modifications were made to the 
bulk data file outside of PATRAN before it was submitted to NASTRAN for analysis. These 
modifications are discussed in this paper. 

Appendices to this work provide additional information about the numeric solutions to 
these problems. Appendix A tabulates NASTRAN’s numeric results for each problem and 
compares it to the appropriate analytic result. MSC/PATRAN was used to produce the fringe 
plots shown. It should be noted that PATRAN uses a linear interpolation between nodes in 
a finite element mesh, no mater what type of element was used. Thus, linear and quadratic 
element meshes having an equal number of nodes will appear identical in a PATRAN-generated 
fringe plot. Appendix B lists the FORTRAN codes written in support of this work. Several 
codes make calls to Bessel function algorithms. These algorithms can be found in the text 
Numerical Recipes: The Art of Scientific Computing (FORTRAN Version fi . Appendix C 
provides computer file usage data for the numerical solutions. 

Section 2: The Cubic Problem 

The fluid/structure interaction problem was first considered for a cubic domain. This 
domain was chosen because it represents a geometrical configuration wherein the wave equation 
is readily solvable in Cartesian coordinates. To further facilitate solution of the equation using 
analytic methods, boundary conditions of P=0 were enforced on all surfaces of the fluid not in 
contact with the structure. For structural portions of the model, thin elastic plates were used 
with boundary conditions of simple-support on each edge. For the forced response analysis, 
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forces of equal magnitude and phase but opposite direction were applied to the model. They 
were applied such that translational vibration modes did not need to be considered for this 
analysis. The combination of these conditions makes an analytic solution of the problem 
straightforward. A skematic representation of the forced response fluid/structure model for 
the cubic geometry is shown in Figure 1. 



Figure 1: Exploded view of the fluid/structure model for the cubic geometry. 

The Free Fluid Problem 

For an ideal, stationary fluid, the acoustic pressure field is described by the wave equation 


V 2 P 


d 2 P _ 
C 2 & 2 " 


( 1 ) 


where, in Cartesian coordinates, 


d 2 P d 2 P d 2 P 


( 2 ) 
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Using the separation of variables technique 2 in conjunction with the homogeneous boundary 
conditions 

P(0,y, z,t) = P(A,y,z,t) = 0 

P(x,0,z,t) = P(x,A,z,t) = 0 (3) 

P(x,y,Q,t) = P(x, y, A,t) = 0 

and assuming a harmonic time dependance, the solution to Equation 1 is 


P(x,y,z,t ) = e iwt X X E D nmk sin—^-sin 

n= 1 m = l k=. 1 


m Try . knz 

-t-™— 


(4) 


where Z) nm * is a constant to be determined from initial conditions. The natural frequencies 
Unmk are given by 


Unmk = \/n 2 + m 2 + Ar 2 (5) 

A cubic fluid volume was discretized with three different meshes for analysis by 
MSC/NASTRAN. One model was constructed using 1000 linear HEX8 elements and 1331 
nodes. The second model used 216 quadratic elements and 1225 nodes. The third model 
was discretized with 1000 quadratic HEX20 elements and 4961 nodes. The fluid in each of 
the models was given the material properties for density and speed of sound shown in Table 
1. Boundary conditions in agreement with Equation 3 were also applied to the each of the 


models. 


Symbol 

Property Name 

Material Property Value 

E 

Young’s Modulus 

10.3 x 10 6 psi 

h 

Thickness of structure 

0.0625 in 

A 

Length of side 

5 in 

C 0 

Acoustic speed of sound 

13.620 x 10 3 in/sec 

V 

Poission’s Ratio 

0.334 

Ps 

Density of structure 

2.5383 x 10~ 4 slugs/in 3 

pi 

Density of fluid 

1.170 x 10 -7 slugs/in 3 


Table 1 Properties for the cubic model. 
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To specify a fluid element within NASTRAN, several bulk data cards must be specified or 
changed. The brief outline of these cards which follows is as shown in the MSC/NASTRAN 
Quick Reference Guide , Version 68. The reader is referred to this source for more detail 
regarding these cards. The first card to be modified is the GRID card. The general format 
of this card is as shown in Figure 2 3 . A “-1” in the CD field of the GRID card is used to 
signify a grid belonging to a fluid element. It should be noted that such a grid point can only 
be defined for volume elements. 


GRID 

ID 

CP 

XI 

X2 

X3 

CD 

PS 

SEID 



ID 

CP 

Xl,X2,X3 

CD 

PS 

SEID 


Grid point identification number. 

Identification number of coordinate system in which the location of the grid point is 
defined. 

Location of the grid point in coordinate system CP. 

Identification number of coordinate system in which the displacements, degrees of 
freedom, constraints, and solution vectors are defined at the grid point. 

Permanent single-point constraints associated with the grid point. 

Superelement identification number. 


Figure 2: NASTRAN GRID bulk data card format. 


The second card that must be specified is the material properties card. The format of the 
material properties card is shown in Figure 3 3 . Note that the MAT10 card specifies material 
properties for fluid elements only. 


MATIO 

MID 

BULK 

RHO 

C 







MID 

Material identification number. 

BULK 

Bulk modulus. 

RHO 

Mass density. 

C 

Speed of sound 


Figure 3: NASTRAN MATIO bulk data card format. 
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The third bulk data card that must be specified is the PSOLID card, shown in Figure 
4 3 . For fluid elements in the model, the FCTN field of the PSOLID card must be specified 
as “PFLUID”. 


PSOLID 

PID 

MID 

CORDM 

IN 

STRESS 

ISOP 

FCTN 




PID Property identification number. 

MID Identification number of a MAT1, MAT4, MAT5, MAT9, or MAT10 entry. 

CORDM Identification number of the material coordinate system. 

IN Integration network. 

STRESS Location selection for stress output. 

ISOP Integration scheme. 

FCTN Fluid element flag. (Character: "PFLUID" indicates a fluid element, "SMECH" 

indicates a structural element; defalut = "SMECH" 

Figure 4: NASTRAN PSOLID bulk data card format. 

MSC/NASTRAN performed a normal modes analysis (SOL 103) on the linear and 
quadratic models. The NASTRAN-calculated frequency for the linear HEX8 model mode 
1,1,1 was 2368.77 Hz. For the 1225-node quadratic HEX20 model, the eigenfrequency for this 
mode was 2359.18 Hz. The 4961-node quadratic HEX20 model yielded a frequency of 2359.32 
Hz. For comparison with NASTRAN’s results, the natural frequencies for these models can 
be determined analytically using Equation 5. From this relation, the first natural frequency 
for each model is 2359.05 Hz, mode 1,1,1. A fringe plot of this mode shape for each of the 
models is shown in Figures Al, A2, and A3. It should be noted that in these (and subsequent) 
modal fringe plots, the displacements have been normalized to one. 

The NASTRAN solution was compared for accuracy against the analytic solution, Equa- 
tion 4. A fringe plot displaying the difference between analytic and numeric solutions at each 
node in the model was created for each case. The error plot for the linear HEX8 cube, mode 
1,1,1 is shown in Figure A4. The error plot for the 1225-node quadratic HEX20 cube, mode 


7 










1,1,1 is shown in Figure A5. Figure A6 shows the error plot for the 4961-node quadratic model. 

This analysis was repeated for mode 1,3,1 of the linear model and mode 1,1,3 of the 
quadratic models. These mode shapes are shown in Figures A7, A8, and A9. The associated 
error plots for each of these mode shapes are shown in Figures A10, All, andAl2. A summary 
of the normal modes analysis for the fluid cube is shown in Table 2. Here, the maximum 
model error is the largest difference anywhere in the model between the normalized analytic 
and numeric solutions. 


Model 

Type 

Elements 

Nodes 

Mode 

Shape 

Eigenvalue 

Error 

Maximum 
Error in 
Model 

Fluid only 

1000 Linear HEX8 

1331 

i,i,i 

0.4117% 

6.83 x 10- 6 

215 Quadratic HEX20 

1225 

0.0055% 

0.0015 

1000 Quadratic HEX20 

4961 

0.0007% 

0.0002 

1000 Linear HEX8 

1331 

1,3,1 

3.134% 

0.510 

215 Quadratic HEX20 

1225 

0.3137% 

0.4593 

1000 Quadratic HEX20 

4961 

0.0432% 

0.386 


Table 2 Results for cubic geometry, normal modes analysis. 


The Free Plate Problem for the Fluid/Structure Model 

Next, the model was modified by bounding two opposing faces of the fluid cube with thin, 
elastic plates. A partially coupled solution to the free vibration problem for this system is 
considered first. For the unforced case, the thin elastic plate obeys the following equation 
of motion 4 : 


DV 4 w + p 3 h 


d 2 w 


= 0 


( 6 ) 


where w is the displacement of the plate in the direction normal to its surface and D, the 
flexural rigidity, is given by 


D = 


Eh 3 

12(l-i/ 2 ) 


(7) 
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Assuming a harmonic time dependence and simply-supported boundary conditions at each 
edge, the out-of plane displacement w can be written 


w(x,y,t) 


e iut 


££«--(“ )-(=?) 
n— 1 m — 1 


with its natural frequencies given by 


^mn — 



( 8 ) 


( 9 ) 


The Free Fluid Problem for the Fluid /Structure Model 

For the ideal fluid between the plates, the wave equation (Equation 1) still applies, but 

with the following boundary conditions: 

P( 0,y,z,t) = P(A,y,z,t) = Q 


P(x,0,z,t) = P(x,A,z,t) = 0 

dP(g,y, ^-,f) _ (Pw(x ,y,<) 

dz P * dt 2 

dP(x,y, f,t) _ d 2 w(x,y,t) 
dz ~ Pf dt 2 

That is, the fluid is constrained to match the behavior of the plates at those points where they 
are in contact. Applying these boundary conditions yields the following expression for P : 


P(x,y,z,t) 


e iwt 


CO OO 

_ . /mnx\ . /niry\ y x 

2^ 2^ D ^ Sm (— J C0S ( a mnZ:) 


where 


_ i 

A 2 ^ 


, „ 2 \ w 
+ n ) ~ rj 


(ii) 


( 12 ) 


and 

r, _ PjUpCmn ('19') 

OmtiMllyttmnjKn ~ w ) 

Three cubic fluid/structure models were constructed for a NASTRAN normal modes 
analysis. The fluid and structure elements within each of the models were given the material 
properties shown in Table 1. The number and type of elements and nodes were as shown in 
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Table 3. The boundary conditions of Equation 10 were applied to the fluid surfaces in the 
model and the edges of the plates were simply supported. 


Model 

Material 

Element Type 

Number of Nodes 

Total Nodes in 
Model 

1 

Structure 

200 Linear QUAD4 

242 

1573 

Fluid 

1000 Linear HEX8 

1331 

2 

Structure 

72 Quadratic QUAD8 

266 

1491 

Fluid 

216 Quadratic HEX20 

1225 

3 

Structure 

200 Quadratic QUAD8 

682 

5643 

Fluid 

1000 Quadratic HEX20 

4961 


Table 3 NASTRAN models for the cubic geometry fluid/structure problem 


For models containing both fluid and structure elements, it is necessary in the NASTRAN 
bulk data file to explicitly specify which nodes lie on the interface between the fluid and 
structure portions of the model. This is accomplished by using the NASTRAN ACMODL 
card. The general format of ACMODL card is shown in Figure 5 3 . For this work, each 
structure grid in the SSET card had a corresponding fluid grid in the FSET card. 


ACMODL 

INTER 

INFOR 

FSET 

SSET 

FSTOL 






INTER Type of interface between the fluid and the structure. 

INFOR Indicates wether FSET and SSET are to be used to define the fluid-structure interface 

FSET Identification number of a SET1 entry that contains a list of fluid grid points on the 

interface. 

SSET Identification number of a SET1 entry that contains a list of structural grid points on the 

interface. 

FSTOL Tolerance, in units of length, used in determining the fluid-structure interface. 

Figure 5: MSC/NASTRAN ACMODL bulk data card format. 

Entry u by hand” of the nodes required for the ACMODL card would be an extraordinarily 
tedious task, particularly for more complicated or larger models. Several FORTRAN codes 
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were created to write this data for the ACMODL card from data already in the bulk data file. 
These codes are listed in Appendix B of this work. 

From Equation 9, the first natural frequency of the plate is 484.54 Hz, mode 1,1. For 
the linear QUAD4 model, the NASTRAN-calculated eigenfrequency for mode 1,1 of the plate 
was 486.86 Hz. For the 266 node quadratic QUAD8 model, it was 481.02 Hz. The 682-node 
QUAD8 model yielded a result of 482.37 Hz. This mode shape is shown for each of the QUAD4 
and QUAD8 models in Figures A13, A14, and A15. From these figures, it is apparent that a 
strong cross-coupling between the front (z=y) and rear (z= — y) plates is occurring. Unlike 
the partially-coupled analytic solution derived above, NASTRAN assumes a fully coupled 
solution. Thus, the first mode shape of the front plate will be coupled through the fluid to the 
rear plate in the NASTRAN solution. In the uncoupled analytic solution, the first mode shape 
of the system is a 1,1 mode shape for one plate and a stationary (0,0) shape for the other. 

The uncoupled analytic solution was used to analyze the mode shapes present on the front 
plate of the model. The difference between the analytic and numeric solutions at each node 
for all three models is shown in Figures A16, A17, and A18. 

This analysis was repeated for plate mode 2,2. This mode shape is shown in Figures 
A19, A20, and A21. The associated error plots are shown in Figures A22, A23, and A24. A 
summary of the normal modes of vibration for the structure portions of the cubic models is 
listed in Table 4. 


Model Type 

Elements 

Nodes 

Mode 

Shape 

Eigenvalue 

Error 

Maximum 
Error in 
Model 


Table 4 Normal modes analysis for structure portion only of the 
fluid/structure model, cubic geometry. (Continued) . . . 
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Fluid/struct 

200 Linear QUAD4 

243 

1,1 

0.4778% 

-1.29 x 10- 3 

72 Quadratic QUAD8 

266 

0.7265% 

-1.82 x 10- 3 

200 Quadratic QUAD8 

682 

0.4481% 

-1.17 x 10~ 3 

200 Linear QUAD4 

243 

2,2 

2.898% 

95.5 x 10“ 3 

72 Quadratic QUAD8 

266 

0.7275% 

133.4 x 10“ 3 

200 Quadratic QUAD8 

682 

0.5628% 

48.9 x 10“ 3 


Table 4 Normal modes analysis for structure portion 
only of the fluid/structure model, cubic geometry. 


For the fluid portions of the system, fringe plots of NASTRAN calculated mode 1,1,0 are 
shown in Figures A25 (linear HEX8 elements), A26 (1225-node quadratic HEX20 model), 
and A27 (4961-node quadratic HEX20 model). For the linear model, NASTRAN calculated 
a frequency for this mode of 1934.09 Hz. For the 1225-node quadratic model, it was 1926.26 
Hz. The 4961-node HEX20 model yielded 1926.17 Hz. A fringe plot showing the difference 
between the analytic and numeric solutions of this mode for each model is shown in Figures 
A28, A29, and A30 

This analysis was repeated for mode 1,1,1 of the fluid for both models. This mode shape 
is shown in Figures A31, A32, and A33. The associated error plots are shown in Figures A34, 
A35, and A36. A summary of the normal modes of vibration for the fluid is listed in Table 5. 


Model Type 

Elements 

Nodes 

Mode 

Shape 

Eigenvalue 

Error 

Maximum 
Error in 
Model 

Fluid/struct 

1000 Linear HEX8 

1331 

1,1,0 

0.4637% 

5.6 x 10~ 6 

216 Quadratic HEX20 

1225 

0.0052% 

728.0 x 10" 6 

1000 Quadratic HEX20 

4961 

0.0007% 

96.1 x 10“ 6 

1000 Linear HEX8 

1331 

1,1,1 

0.4117% 

6.6 x 10" 6 

216 Quadratic HEX20 

1225 

0.0054% 

1.532 x 10- 3 

1000 Quadratic HEX20 

4961 

0.0007% 

198.6 x 10- 6 


Table 5 Normal modes analysis for fluid portion 
only of fluid/structure model, cubic geometry. 
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The Forced Fluid /Structure Problem 

The third and final analysis performed for the cubic geometry model was that of forced 
vibration. A finite element model identical to the linear model used for the free vibration 
analysis was used, the only difference being the application of a harmonic point force at the 
center of each plate in the system. The directional sense of these forces was such that they 
pointed in toward the fluid volume between the plates. Again, the analytic solution of the 
uncoupled problem is considered first. For a forced vibration analysis of a thin, elastic plate 
including damping effects, the governing equation is now written 

DV 4 w + 7 ^ + pillar = F(x,y,t) (14) 

where 7 is the viscous damping coefficient and F( x,y,t) is a harmonic point force applied to the 
plate. For the case where this force is applied to the center of the plate, it can be written as: 

F(x,yJ) = F 0 e iu ' , 6 (x-iy(y-^j (15) 

Applying simply-supported boundary conditions to the plate, the solution to Equation 14 is 
given by 

oo oo 

£ £c mn »m(!^).m(^) ( 16 ) 

m= 1 n = l 

Expanding F(x,y,t) in terms of the eigenfunctions of the plate, we have 

OO CO 

F( I , > ,l) = a'-'££f m ^m(^)»in(^) (17) 

m—0 n=0 

where 

F mn = ^- 08) 

and thus 


4 F P sm(^)sm(qf) 

A 2 p a h -Ul 2 + 2j»/U>U) mn ) 


(19) 
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Note that in Equation 16, the viscous damping coefficient y has been replaced by a frequency- 
dependent damping coefficient tj. The relation between these two coefficients is: 


p,h ~ 2T7Wmn ( 20 ) 

For the fluid in the region between the plates, Equation 1 and boundary conditions 10 still 
apply. Using the above results for the plates, the equation for acoustic pressure at a point 
in the fluid is now written 


CO oo 

P(x,y,z,t) = -e' u, ^T D mn sin(J^ S jsin(J^^cos(a mn z) 


where again 


9 7T / o 9\ UJ 

a mn ~ 72 ( m +n ) ~ Z2 (22) 

c o 

but 

D _ 4F„ 

A2 P> h Qmr l sm(^ 1 )(w2 in - u) 2 + 2irfjiui mn ) ^ 

For the numerical analysis by NASTRAN of this forced- response problem point forces 
of -5sin(u>t)n lbs. were applied to the center of each plate in the linear model. Boundary 
conditions for the fluid and structure portions of the model were otherwise identical to those 
used for the free vibration problem. 

A direct frequency response analysis (SOL 108) was performed by NASTRAN over a 
frequency range from zero to 5000 Hz. Structural damping in the model was specified by using 
the PARAM G card in the NASTRAN bulk data file. Figure A37 shows the displacement of a 
node located at the center of the front plate for the linear QUAD4 model. Figure A38 shows 
the acoustic pressure disturbance at a node halfway between the center of the front plate and 
the center of the fluid region. Also shown are the partialy-coupled analytic solutions for the 
displacement and pressure respectively at these two locations. 
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A modal frequency response analysis (SOL 111) was next performed by NASTRAN on 
the same cubic-geometry model. The first twenty NASTRAN-calculated fluid and structural 
modes were used by NASTRAN to calculate the overall response of the model. It should 
be noted that the twentieth natural frequency of the structure, as calculated by NASTRAN, 
occurs at approximately 3000 HZ, while the twentieth natural frequency of the fluid occurs 
above 5000 Hz. The results of this analysis for a frequency range of zero to 5000 Hz are shown 
in Figures A39 for a point at the center of the front plate, and A40 for a point midway between 
the front plate and the center of the cube, along with analytic solutions at each of these points. 

Section 3: The Cylindrical Problem 

A cylindrical domain was next defined for analysis. As was the case for the cubic domain, 
a cylindrical geometry was chosen because an analytic solution of the wave equation in 
cylindrical coordinates becomes straightforward and uncomplicated. To facilitate solution 
of the wave equation using the separation of variables technique, boundary conditions of P=0 
were again enforced on all fluid surfaces not in contact with any structures in the model. 
For models in which a structure was present, a thin cylindrical shell was used. Boundary 
conditions for the shell were chosen which simplified the analytic solution. For the forced 
response analysis, forces applied to the model were equal in magnitude and phase, but opposite 
in direction. They were applied to the cylindrical shell such that translational vibration modes 
would not have to be considered in this analysis. A skematic representation of the cylindrical 
fluid/structure geometry for the forced response analysis is shown in Figure 6. 
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Fluid/Structure Cylindrical Geometry, Exploded View 



Figure 6: Exploded view of fiuid/structure model for the cylindrical geometry. 


The Free Fluid Problem 


For an ideal, stationary fluid, the wave equation in cylindrical coordinates now becomes 

<*> 

where V 2 P is given by 


d 2 P 1 dP 1 d 7 P d 2 P 


V^ = t-t + -w- + 


+ 


dr 2 r dr r 2 dd 2 dz 2 

Assuming a harmonic time dependance and homogenous boundary conditions: 

P(r,$,0) = P(r,e,l) = 0 


(25) 


P(a,B,z)= 0 

|P(O,0,z)| < oo (boundedness) 

P(r, 6, z) = P(r, 9 + 27rn, z) (periodicity) 
separation of variables yields the solution 

OO OO OO 

P(r,9,z,t) = e lwt ^ Jr l (J^^sin(J-^^^[B mnj cos(n8) + C mn jsin(n9)] 

j= 1 m — 1 n = 0 


(26) 


(27) 
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where J n is the nth order Bessel function, rj is its jth root, and a and l are the radius and 
length of the cylinder, respectively. The natural frequencies uj m n are given by 

Three cylindrical fluid volumes were discretized for analysis by MSC/NASTRAN. One 
model was constructed using 2240 linear WEDGE6 elements and 1449 nodes, one used 348 
quadratic WEDGE15 elements and 1200 nodes, and the third used 2240 quadratic WEDGE15 
elements and 6609 nodes. The fluid in each of the models was given the material properties 
for density and speed of sound shown in Table 6. Geometric dimensions are also shown in this 
table. Boundary conditions in agreement with Equation 26 were applied to each model. 


Symbol 

Property Name 

Material Property Value 

E 

Young’s Modulus 

10.3 x 10 6 psi 

a 

Radius of cylinder 

1 in 

Co 

Acoustic speed of sound 

13.620 x 10 3 in/sec 

h 

Thickness of structure 

0.0625 in 

l 

Length of the cylinder 

5 in 

V 

Poission’s Ratio 

0.334 

Ps 

Density of structure 

2.5383 x 10 -4 slugs/in 3 

Pi 

Density of fluid 

1.170 x 10” 7 slugs/in 3 


Table 6 Properties for cylindrical model. 


For each model, NASTRAN performed a normal modes analysis (SOL 103). The natural 
frequencies for this geometry can be determined analytically from Equation 28 for comparison 
to NASTRAN’s results. Analytically, the first natural frequency for this model is 5387.91 
Hz, mode 1,0,1. The NASTRAN linear WEDGE6 model calculated an eigenfrequency of 
5445.67 Hz for this mode. The quadratic 348-element model yielded a frequency 5392.38 Hz. 
The mode 1,0,1 eigen frequency for the 2240-element WEDGE15 model was 5388.06 Hz. A 
normalized fringe plot of this mode shape is shown for each model in Figures A41, A42, and 
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A43. The NASTRAN numeric solution was compared to the analytic solution, Equation 27. 
A fringe plot for each model showing the deference between the analytic and numeric solutions 
at each node is shown in Figures A44, A45, and A46. 

This analysis was repeated for mode 1,1,3 for each model. This mode shape is shown in 
Figures A47, A48, and A49. However, at the time of this writing, it was not possible to create 
an error fringe plot for this mode. In calculating the eigenvectors for a model in which a 
rotational symmetry is present, NASTRAN introduces an arbitrary phase angle with respect 
to the analytic solution into its numeric solution. That is, the numeric solution is accurate, 
but its mode shape is rotated by some angle about its axis of symmetry. It was not possible to 
adequately determine what this angle was. Thus, a node-by node comparison to the analytic 
solution was not possible. A summary of the cylindrical normal modes analysis is shown in 
Table 7. The maximum model error is again the largest difference anywhere in the model 
between the normalized analytic and numeric solutions. 


Model Type 

Elements 

Nodes 

Mode 

Shape 

Eigenvalue 

Error 

Maximum 
Error in 
Model 


2240 Linear WEDGE6 

1449 


1.07% 

0.017 


348 Quadratic WEDGE15 

1200 

1,0,1 

82.87 x 10" 3 % 

0.005 

Fluid only 

2240 Quadratic WEDGE15 

6609 


6.50 x 10~ 3 % 

774.2 x 10 -6 

2240 Linear WEDGE6 

1149 


3.04% 

N/A 


348 Quadratic WEDGE15 

1200 

1,1,3 

0.356% 

N/A 


2240 Quadratic WEDGE15 

6609 


0.015% 

N/A 


Table 7 Results for cylindrical geometry, normal modes analysis. 


The Free Shell Problem 


The next problem considered for the cylindrical geometry was that of a fluid-structure 
interaction. We begin by defining a thin, finitely-long cylindrical elastic shell filled with a 
stationary, ideal fluid. As we did for the cubic fluid/structure geometry, we shall solve the free 
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vibration problem for this system using an uncoupled solution. For a thin, cylindrical shell, 
the Donnell-Mushtari equations of motion are: 5 


where 


cl) di 2 


/ 1 \ d 2 v / 1 + tA ' 
\ 2a )l 


d 2 u j 

M + i/> 

\ d2 U f 

i + d 2 v 

v dw 

— n 

" \ 

< 2a 2 ) 

1 d9 2 V 

2 a 2 ) dzdO 

a dz 

— u 

d 7 u 

(\-v' 


1 \ d 2 v ( 

1 \ dw 

— n 

dzde 1 

K 2 , 

) dz 2 v 

a 2 ) d0‘> V 

)~de 

— u 

\d 2 tc 
) dt 2 + 

v / du \ / 1 \ 

a\dz) \a 2 / 

dv 1 

^ H — ~ w + 
dv cr 


= 0 


(29) 


C L = 


M 


(30) 


On the LHS of Equation 29, u(z, 6), v(z,8), w(z,9) represent displacements in the axial, 
circumferential, and radial directions respectively. Assuming a harmonic time dependence 
with the following boundary conditions 5 : 

«(0, 0) = v(l,0) = 0 

(31) 


u>(0, 0) — w(l, 9) = 0 


die complete solution to Equation 29 can be written as: 


u(9, z ) = e iu,t Y) Y2 C0S (^J~) [A mn cos(n0) + ^„sin(n«)] 

ri=rO m = 0 

oo oo 

v{6,z) = + (32) 

n=0 n= 0 

oo oo 

w(0, z) = e iwt Y2 sin (^j- L )[C mn cos(n0) + C^„sin(nfl)] 

n=0 m — 0 

Because of the boundary conditions imposed for this analysis, either the starred or un-starred 
terms can be considered as a complete solution 6 . This work will use the un-starred solution. 
That is, we will take 


OO OO 

u{z,6) = e‘"‘ Y2 Y2 A mn cos ^ ^ cos(n6) 

n=0 m=0 

oo oo 

v(zj) = e iwt ^2 J2 B rnnSin\J^jsin(nO) (33) 

n — 0 m = 0 

oo oo 

w(z, 9) = e iwi Cmn s * n ( —y~) cos(nO ) 
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The end conditions shown in Equation 31 represent “shear diaphragms" 7 . Non — axial dis- 
placements at the ends of the cylindrical shell (at z = 0 and z = /) are zero. Substituting 
Equation 33 into Equation 29 yields: 


A-f^An 


(-A 2 - i^An 2 + QlJ 

iiiiAn (_i^A 2 -n 2 + ^J 


t/A 

— n 


— v\ 


(l + T^(> 2 + n 2 ) 2 -^„)JLC, 



^mn 


o' 


Bmn 

= 

0 


Cmn 


0_ 


where 


and 


2 _P>{ l-^ 2 )a 2 w 2 


= 


A = 


mna 


(34) 


(35) 


(36) 


For non-trivial solutions to Equation 34, the determinant of the coefficient matrix is set equal 
to zero. Doing so produces the following characteristic equation in fl^j n 7 : 


^ mn -K^ mn +K^ mn -K G = Q 


(37) 


where 


K 2 - 1 + —^—(n 2 + A 2 ) + —j (n 2 + A 2 ) 2 


A'i - — — (3 + 2t/)A 2 + n 2 + (n 2 + A 2 ) + ^ (” 2 + ^ 2 ) 3 


(38) 


A'n = 


1 - v 


(l-i/ 2 )A< + 


h 2 

12a 2 


(n 2 + A 2 / 


The roots of this equation in combination with Equation 35 can be used to calculate the 
natural frequencies uj mn of the cylindrical shell. 

A more accurate description of the motion of thin cylindrical shell is provided by the 
Epstein-Kennard equations of motion. Unlike the Donnell-Mushtari formulation, the Epstein- 
Kennard relation does not neglect changes in shear stresses acting through the thickness of 
the shell. For a more detailed analysis of this particular shell theory, the reader is referred to 
the literature [7]. Here, it is noted that, with respect to the natural frequencies of vibration, 
Equation 37 becomes 7 


- ( A 2 + kAI< 2 )Q 4 mn + (A’i + kAk\)n 2 mn - (A'o + kAKo) = 0 


(39) 
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where, for the Epstein-Kennard theory, 

1-1-31/ (2 - 8i/ 2 + 3 i/ 3 )A 2 (19 - 37i/ + 19i/ 2 + i/ 3 ) i/ 2 (n 2 + A 2 ) 

AAz= l-v 2(1 - vf 2(1 — i/) 2 (l-*) 2 

(3 + 8u - 5i/ 2 - i/ 3 )A 2 (2 + i/)n 2 (6 + 4i/ - 8i/ 2 + 3i/ 3 - 8i/ 4 )A 4 i/ 2 (n 2 + A 2 )" 

AKl= 2(l-i/) + 2 4(l-i/) 2(1-!/) 

(26- 60i/ + 40i/ 2 -3i/ 3 - 8i/ 4 )A 2 n 2 (l3-22i/+ 10i/ 2 )n 4 

2(1 -!/) 2(1-j/) 

(2 + 6t/ - 2t/ 2 - 3i/ 3 ) A 4 | 4A2n2 f n4 1 + «/ >6 7-5v x4 _ a cx2 „ 4 


( 40 ) 


AK 0 = 5(1-1/) 


2 ( 1 -!/) 


1 -1/ 


-A” — 


1 - 1/ 


-A 4 n 2 -8A 2 n 4 -2n* 


and 


k = 


/i 2 

12a 2 


(41) 


and the natural frequencies can be calculated as before. 

The Free Fluid Problem 

For the fluid-filled region inside the cylindrical shell, the wave equation in cylindrical 
coordinates (Equation 24) can be used. The boundary conditions now become 

P(r,0,O) = PM,/) = O 


dP(a,6,z) d 2 w(0,z) 

d r = Pl W- 

|P(O,0,r)| < 00 (boundedness) 


(42) 


P(r, 6, z) = P(r, 0 + 2 nn, z) (periodicity) 

where u;(0, c) represents the radial displacement of the cylindrical shell. That is, just as in the 
case for the cubic geometry, the fluid is constrained to match the behavior of the structure 
at those points where they are in contact. Applying these boundary conditions and using the 
results for the cylindrical shell yields the following expression for P 5 : 


OO OO 

P{r,0 , z) = e'"‘ ^2 DmnJn{ a mn r ) s in(J—j—}COs(n8) 


(43) 


n— 0 m = l 


where 


Pm n — 


C m nP/W 2 


Q mn ^‘*Ai( t *mn a ) ^mn^n + l(^mn a )j 


(44) 
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and 


2 (Umn \ 2 ( mn \ 2 

—) -(—) (45) 
Three cylindrical fluid/structure models were constructed for NASTRAN normal modes 
analysis. The three models were discretized as shown in Table 8. The fluid and structure 
elements within each of the models were given the material properties shown in Table 6. 
Boundary conditions as shown in Equations 31 and 42 were applied on each model. 


Model 

Material 

Element Type 

Number of 
Nodes 

Total Nodes in 
Model 

1 

Structure 

480 Linear QUAD4 

504 


Fluid 

2240 Linear WEDGE6 

1449 

1953 

9 

Structural 

192 Quadratic QUAD8 

608 

2729 

£ 

Fluid 

672 Quadratic WEDGE15 

2121 

Q 

Structure 

480 Quadratic QUAD8 

1488 

8097 

o 

Fluid 

2240 Quadratic WEDGE15 

6609 


Table 8 NASTRAN models for fluid/structure problem, cylindrical geometry. 


From Equation 35, using Epstein-Kennard theory, the natural frequency of the cylindrical 
shell for mode 1,1 is 6248.99 Hz. Using the linear QUAD4 model, NASTRAN calculated an 
eigenfrequency of 6251.14 Hz. The 192-element QUAD8 model yielded a frequency for this 
mode of 6256.34 Hz. A frequency of 6245.55 Hz was calculated using the 480-element QUAD8 
model. This mode shape is shown for the u, v, and w displacements for each model in Figures 
A50, A51, A52, A53, A54, A55, A56, A57, and A58. The “breathing modes” of the shell, 
i.e. modes where m=l and n=0 in Equation 33 occur (using Epstein-Kennard theory) for this 
model at 12,318.19 Hz, 19,585.91 Hz and 35,037.39 Hz. The second frequency, 19,585.91 Hz, 
is shown for each model in Figures A59, A60, A61, A62, A63, A64, A65, A66, and A67, with 
error plots for the radial displacements in Figures A68, A69, and A70. 
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For the fluid inside the cylindrical shell, fringe plots of NASTRAN-calculated mode 1,0,1 
is shown in Figures ATI (linear WEDGE6 elements), A72 (672-element WEDGE15 model) 
and A73 (2240 element WEDGE15 model). A fringe plot showing the difference between 
the analytic and numeric solutions of this mode is shown in Figures A74, A75, and A76. 
A summary of the normal modes analysis for the structure portion of the fluid/structure 
cylinder is shown in Table 9. 


Model Type 

Elements 

Nodes 

Mode 

Shape 

Eigenvalue 

Error 

Maximum 
Error in 
Model 


2240 Linear WEDGE6 

1449 


0.102% 

98.4 x 10~ 6 


672 Quadratic WEDGE15 


1,0,1 

2.937 x 10~ 4 % 

0.0824 


2240 Quadratic 
WEDGE15 



101.3 x 10“ 6 

Fluid/struct. 

480 Linear QUAD4 

504 


0.143% 

N/A 

192 Quadratic QUAD8 


1,1 

0.118% 

N/A 


480 Quadratic QUAD8 

1488 


0.055% 

N/A 


480 Linear QUAD4 

504 


0.397% 

0.087 x 10~ 6 


192 Quadratic QUAD8 

608 

1,0 

0.294% 

-0.0270 


480 Quadratic QUAD8 

1488 


0.040% 

13.279 x 10" 3 


Table 9 Results for cylindrical fluid/structure geometry, normal modes analysis. 

The Forced Fluid /Structure Problem 

The third problem considered for the cylindrical fluid/structure model was that of a forced 
response. As we did for the cubic model, we first consider the analytic partially-coupled 
solution of this problem. The cylindrical shell equations of motion with viscous damping are 
given by 5 : 

L n L12 
L>21 L22 

7*31 ^32 


1 

cn 


U 

L>23 


V 

^33. 

DM 

IV 


l-I / 2 

Eh 


fz ' 

fe 

L"/r. 


(46) 
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where 


L 11 
^12 
^13 
^21 
£ 22 
^23 
^31 
^32 

£33 


= -w 2 + 2i 

_ _ 1 + i' 3 2 

2a 3z30 

_ ^3 

a 3.2 

1 + v 3 2 

2 a 3230 

1 a 2 (I-*/ 2 ) 5 1 -* a 2 1 a 2 

C 2 L dt 2+ ' r Eh m 2 d^~~E*W 

= -—— 
a°~de 

_ V d 
a dz 

= -LA 

a 2 30 

1 a 2 ( 1 -* 2 ) a 1 h 2 , 

- pry a 72 + 7 — 777 — -77 + — 4 - 777 V 4 

G£ ot 2 £7i a< a 2 12 


( 47 ) 


and /*> /»> and / r represent forces in the axial, circumferential, and radial directions. Here 
the convention of a positive inward-directed radial force has been assumed^. Applying forces 
to the cylinder of 


fz =0 


fe = 0 


f r = F 0 sin(u>t)6 



5) [*(*) + *(«-»)], 


( 48 ) 


expanding these forces in terms of the eigenfunctions of the shell and once again using Equation 
33 with the boundary conditions of Equation 31, Equation 46 yields 


Ln L\2 ^13 


A-mri 


' 0 

Z/21 L 22 ^23 


Fmn 

1 

A L 

0 

_^31 £32 £33 _ 


J-'mn _ 

Psfl 

fmn_ 


(49) 
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where now 


1 - u 


L\\ — — u + 2 irjujuj mri - 1 - -f — — n 2 Cj v 


L 12 — 2 1 — — 




^13 — ^31 — 


L 22 — + 2.irju)u) rnn H — — ^C\ + ~^2^l 


L23 = £32 — —Cl 
0 “ 

L 33 = -w 2 4 2 ir/woJ rnfl + ^ + n 2 ) 2 


A = 


m 7r 


( 50 ) 


(51) 


and F mn is given by 


4 . { mn 

F mn = -F 0 s t n{ — 


) [cos(nTr) + 1 ] 


(52) 


Note that in these equations, the viscous damping term 7 has once again been replaced by 
a frequency-dependent damping term 77. 

For the fluid inside the shell, we proceed as we did for the undamped, free vibration case. 
That is, applying the boundary conditions of Equation 42 to the cylindrical wave equation 
(Equation 24), the solution for the acoustic pressure field becomes : 5 

P(r,e,z) = e iut Y D mn J n (a m „r)sin(—j-)cos(n 6 ) (53) 


with 


and 



(54) 


(55) 


For the NASTRAN analysis of this problem, a cylindrical finite element model was 
constructed using linear elements. The finite element mesh was identical to the linear mesh 



used for the cylindrical free vibrations problems. The geometric characteristic of this problem 
were modified as shown in Table 10. The structural and fluid elements were given the material 
properties shown in Table 6. 


Radius (a) 

10.0 in 

Length (1) 

50.0 in 

Shell Thickness (h) 

0.0625 in 


Table 10 Geometric dimensions for the cylindrical 
geometry forced-response analysis model. 

Boundary conditions in agreement with Equations 31 and 42 were applied to the appro- 
priate structural and fluid elements respectively. Two point forces of 50 sin ( u>t)k lbs., in phase 
and both pointing radially inward were located at 6 = 0, z = | and 6 = ir, z = NASTRAN 
performed a direct frequency response (SOL 108) over the frequency range zero to 5000 Hz. 

Analytic and NASTRAN-generated numeric solutions for the frequency range zero to 1000 
Hz are shown in Figures A78 and A77. Figure A78 shows the displacement vs. frequency of 
a structural node located at r = a, 6 = 0 , z = j. Figure A77 shows the acoustic pressure of 
a fluid node at location r = |, 6 = 0, z = ^. 

Section 4: Conclusions 

MSC/NASTRAN has been used to numerically calculate a number of different mode 
shapes for a cubic and a cylindrical acoustic cavity. The accuracy of the NASTRAN fully- 
coupled fluid/structure solution improves with the use of elements having a greater number of 
nodes (i.e. linear vs. quadratic elements). The difference in performance between these two 
elements becomes particularly significant in models involving a curved geometry. This result is 
not unexpected, as, it is difficult to model accurately a curved shape using only straight linear 
elements. No substantial change in performance was noted in going from a fluid-only model 
to a model containing both fluid and structure elements. NASTRAN was able to compute the 


26 






normal modes for each model quickly and accurately. In sum, NASTRAN’s calculation of the 
normal modes (SOL 103) for a model tended to be simple, direct, quick, and accurate. 

Hard disk space limitations were factor that became an important consideration during 
the forced response analysises, particularly when the direct method (SOL 108) was used. In 
general, rather than use system memory, NASTRAN writes data to files during the solution 
of a finite element problem. Although most of these files are deleted when the NASTRAN 
solution is completed, they can become quite large during the process of solving the problem. 
If enough free space is not available for use by NASTRAN, the problem cannot be solved. 
To work around this problem, the forced response analysis for the cylindrical geometry was 
submitted to NASTRAN as four separate problems, each covering a range of 250 frequencies. 
After the solutions for each problem were calculated, they were combined to form a complete 
solution for the entire frequency range of interest. Appendix C of this work contains a listing 
of the significant files produced by NASTRAN during each analysis and the size of each file. 
While this method was workable for the problem at hand, it could become very unwieldy with 
the addition of more nodes and/or elements to the model. 

In short, NASTRAN’s normal modes analysis was very robust. For a reasonable model 
discretization, its results can be considered accurate. However, for the cases where a forced 
response analysis of a fluid/structure model is required, careful consideration should be given 
to ways in which the model can be simplified and which frequencies are of prime interest. 
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Appendix A Figures 
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Figure Al: Mode 1,1,1 for cubic geometry; 1000 linear HEX8 elements, 1331 nodes. 



Figure A2: Mode 1,1,1 for cubic geometry, 
215 quadratic HEX20 elements, 1225 nodes. 
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Figure A5: Mode 1,1,1 error for cubic geometry, 
215 quadratic HEX2Q elements, 1225 nodes. 
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Figure A6: Mode 1,1,1 error for cubic geometry, 
1000 quadratic HEX20 elements, 4961 nodes. 
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Figure A9: Mode 1,1,3 for cubic geometry, 
1000 quadratic HEX20 elements, 4961 nodes. 



Figure A 10: Mode 1,3,1 error for cubic geometry, 
1000 linear HEX8 elements, 1331 nodes. 
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Figure A 15: Mode 1,1 for structure portion of fluid/structure 
cube. 200 quadratic QUAD8 elements, 682 nodes. 
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Figure A16: Mode 1,1 error for cubic fluid/structure geometry (front 
plate only). 200 quadratic QUAD4 elements, 242 nodes. 
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Figure A 17: Mode 1,1 error for cubic fluid/structure geometry (front 
plate only). 72 quadratic QUAD8 elements, 266 nodes. 
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Figure A18: Mode 1,1 error for cubic fluid/structure geometry 
(structure only). 200 quadratic QUAD8 elements, 682 nodes. 
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Figure A21: Mode 2,2 for structure portion of 
fluid/structure cube. 200 QUADS elements, 682 nodes. 
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Figure A22: Mode 2,2 error for cubic fl u i d/s truct lire geometry 
(front plate only). 200 linear QUAD4 elements, 242 nodes. 
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Figure A24: Mode 2,2 error for cubic fluid/structure geometry (front 
plate only). 200 quadratic QUADS elements, 682 nodes. 
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Figure A25: Mode 1,1,0 for fluid portion of fluid/structure 
cube. 1000 linear HEX8 elements, 1331 nodes. 
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Figure A26: Mode 1,1,0 for fluid portion of fluid/structure 
cube. 216 quadratic HEX20 elements, 1225 nodes. 
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Figure A27: Mode 1,1,0 for fluid portion of fluid/structure 
cube. 1000 quadratic HEX20 elements, 4962 nodes. 



Figure A28: Mode 1,1,0 error for fluid/structure geometry 
(fluid only). 1000 linear HEX8 elements, 1331 nodes. 
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Figure A29: Mode 1,1,0 error for fluid/structure geometry 
(fluid only). 216 quadratic HEX20 elements, 1225 nodes. 
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Figure A30: Mode 1,1,0 error for cubic fluid/structure geometry 
(fluid only). 1000 quadratic HEX20 elements, 4962 nodes. 
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Figure A33: Mode 1,1,1 for fluid portion of fluid/structure 
cube. 1000 quadratic HEX20 elements, 4962 nodes. 
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Figure A34: Mode 1,1,1 error for cubic fluid/structure 
geometry (fluid only). 1000 linear HEX8 elements, 1331 nodes. 
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Figure A36: Mode 1,1,1 error for cubic fluid/structure geometry 
(fluid only). 1000 quadratic HEX20 elements, 4962 nodes. 
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W*Displacement vs. Frequency 

Flukl/structure cube. 1200 Elements 1573 Nodes, NASTRAN direct frequency response analysis 



Figure A37: Displacement at the center of the z=5 plate 
on the fluid/structure cube. NASTRAN direct frequency 
response analysis. Analytic and numeric solutions shown. 


Acoustic Pressure vs. Frequency 
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Figure A38: Acoustic pressure at the point (2. 5, 2. 5, 3. 5) 
in the fluid/structure cube. NASTRAN direct frequency 
response analysis. Analytic and numeric solutions shown. 
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W- Displacement vs. Frequency 

Fkjid/structura cU>* 1200 Elemanta, 1573 Nodas. NASTRAN modal traquancy raaponaa aruOyaia 



Figure A39: Displacement at the center of the z=5 plate on 
the fluid/structure cube. NASTRAN modal frequency 
response analysis. Analytic and numeric solutions shown. 
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Figure A40: Acoustic pressure at the point (2. 5, 2. 5, 3.5) in 
the fluid/structure cube. NASTRAN modal frequency 
response analysis. Analytic and numeric solutions shown. 
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Figure A43: Fluid mode 1,0,1 for cylindrical geometry, 
2240 quadratic WEDGE15 elements, 6609 nodes. 



Figure A44: Fluid mode 1,0,1 error for cylindrical 
geometry, 2240 linear VVEDGE6 elements, 1449 nodes. 
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Figure A45: Fluid mode 1,0,1 error for cylindrical 
geometry, 348 quadratic WEDGE 15 elements, 1200 nodes. 
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Figure A46: Fluid mode 1,0,1 error for cylindrical 
geometry, 2440 quadratic WEDGE15 elements, 6609 nodes. 
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Figure A49: Fluid mode 1,1,3 for cylindrical geometry, 
2240 quadratic WEDGE15 elements, 6609 nodes. 



Figure A50: U-displacement (axial) of cylindrical shell, 
mode 1,1. 480 linear QUAD4 elements, 504 nodes. 
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*- 1.000 * 

Figure A5i: V-displacement (circumferential) of cylindrical 
shell, mode i,i. 480 linear QUAD4 elements, 504 nodes. 



Figure A52: VV-displacement (radial) of cylindrical shell, 
mode 1,1. 480 linear QUAD4 elements, 504 nodes. 
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Figure A55: YV-displacement (radial) of cylindrical shell, 
mode 1,1. 192 quadratic QUADS elements, 608 nodes. 
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Figure A56: U-displacement (axial) of cylindrical shell, 
mode 1,1. 480 quadratic QUAD8 elements, 1488 nodes. 
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Figure A59: U-displacement (axial) of cylindrical shell, mode 1,0 
(breathing mode). 480 linear QUAD4 elements, 504 nodes 
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Figure AGO: V-displacement (circumferential) of cylindrical 
shell, mode 1,0 (breathing mode). 480 linear QUAD4 
elements, 504 nodes. Eigenvector not normalized to one. 
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Figure A63: V-displacement (circumferential) of cylindrical shell, mode 1,0 
(breathing mode). 192 quadratic QUAD8 elements, 608 nodes. 



Figure A64: YV-displacement (axial) of cylindrical shell, mode 1,0 
(breathing mode). 192 quadratic QUADS elements, 608 nodes. 
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Figure A67: YV-displacement (radial) of cylindrical shell, mode 1,0 
(breathing mode). 480 quadratic QUADS elements, 1488 nodes. 
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A68: YV-displacement (radial) error for cylindrical shell, 
(breathing mode). 480 linear QUAD4 elements, 504 nodes. 
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Figure A70: YV-displacement (radial) error for cylindrical shell, mode 1,0 
(breathing mode). 480 quadratic QUAD8 elements, 1488 nodes. 
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Figure A72: Acoustic pressure inside cylindrical shell, fluid mode 
1,0,1. 672 quadratic VVEDGE15 elements, 2121 nodes. 





(fluid only). 2241 linear YVEDGE6 elements 1449 nodes. 
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Figure A76: Mode 1,0,1 error for cylindrical fluid/structure geometry (fluid 
portion only). 2241 quadratic WEDGE 15 elements, 6609 nodes. 
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Numeric and Analytic Displacements vs. Frequency, Cylindrical Geometry 

Radial Displacement at location (x,y,z) a,0,L/2 (node 283) 



Frequency (Hz) 


Figure A77: Radial displacement at the coordinates r = a, 0=0, z = - for 
the fluid/structure cylinder. NASTRAN direct frequency 
response analysis. Analytic and numeric solutions shown. 
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Acoustic Pressure (log) 


Numeric and Analytic Acoustic Pressure Fields vs. Frequency, Cylindrical Geometry 

Acoustic Pressure at location (x t y,z) a/2,0,L/2 {node 1218) 



Figure A78: Acoustic pressure at the coordinates 
r = 9 — 0, z = £ for the fluid/structure cylinder. NASTRAN 

direct frequency response. Analytic and numeric solutions shown 



Appendix B: FORTRAN Codes 


This Appendix contains a source code listing of the FORTRAN programs written in 
support of this work. A brief description of each of these codes is shown in Table 11. 


General Data Manipulation 

prefluid .f 

NASTRAN bulk data file is modified such that it contains fluid elements 
only. 

pch2res.f 

Converts NASTRAN .pch file to a PATRAN-readable .res file. Fluid 
pressures at each node are written as "displacements" in the x-direction. 

convert .f 

Converts NASTRAN .pch file to a TECPLOT data file. For use with 
forced-response analysis 

Written for Cubic Geometry Models 

reader. f 

Compares analytic and numeric normal modes solutions for fluid-only cube 
at each node. 

setmaker.f 

Creates ACMODL card for a cubic geometry. 

plate.f 

Compares analytic and numeric normal modes solutions for the 
fluid/structure cube at each node. 

forres.f 

Computes analytic forced-response at a given location for the fluid/structure 
cubic geometry. 

Written for Cylindrical Geometry Models 

compare.f 

Compares analytic and numeric normal modes solutions at each node for 
fluid-only cylinder. 

card.f 

Creates ACMODL card for a cylindrical geometry. 

modesm.f 

Calculates the natural frequencies of a thin, elastic, cylindrical shell using 
Epstien-Kennard theory. 

forced2.f 

Computes analytic forced-response at a given location for the fluid/structure 
cylindrical geometry. 


Table 11 Description of FORTRAN codes used. 


The codes listed in this Appendix have been written specifically for this study and should 
be considered to be research code only. They are provided for completeness. Their successful 
operation cannot be guaranteed outside the scope of this work. 
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Program prefluid. f 

C1234412 3456789012 

program pre_fluld 

character*8 minusl, fstB, snd8, trd8, frt8, fth8, six8, 

& Sth8,eth8,nth8, tth8 
character* 20 zfile.ofile 
minusl* ' -1 

write(*„ *) 'ENTER BDF FILE NAME ' 
r«ad ( * , 1 ) zfile 
1 format (a) 

of ile= ' P' //zfila 

op#n(unit*10, file=zf ile, status* ' old' ) 
opan (unit=ll , fil#=ofile, status= 'unknown' ) 

100 raad { 10 , 1121 , and=200 ) f st8 , snd8 , trd8 , f rt8 , fth8, six8 
4 , sth8, eth8,nth8, ttta8 

if (fstS.eq. 'GRID '(than 

writ® (11 , 1121 ) f st8 , snd8 , trd8 , f rtB , f th8 , six8 , minusl , 

4 ath8 , nth8, tth8 

alaaif (fst8.eq. 'MAT1 ' ) than 
fstSs'MATIO 
trd8= ’ 

frt8=’ 1. 170E-7' 
fth8=' 13620.0 ' 

writ# (11, 1121) £st8, snd8, trd8, £rt8, fth8, six8 
alaaif (fstS.aq. 'PSOLID ' ) than 
frt8= ' 

£th8= ' 
six8= ' 
sth8= ' 

eth8='PFLUID ' 

nth8= ' 

tth8s' 

writ* (11, 1121) fst8, snd8, trd8, frt8, fth8,six8, sth8, 

4 athS , nthS, tth8 

alsaif (fst8(l: 6) .aq. 'ASSIGN' ) then 
alsaif ( f st8 (1:4) .aq. 'TIME' ) than 
writadl, *) 'TIME 600’ 
do 110 ii*l, 600 

raad (10. 1121, end=200) fst8, snd8, trd8, frt8, fth8, sixS 
4 , sth8 , athS , nth8 , tth8 

if (fstS (1:4) . aq. 'CEND' ) goto 120 
110 continue 

120 writadl, •> 'CEND' 

alsaif ( fstS.aq. ' KETHO'Jthan 
writadl, * } f st8 . 'D(STRUCT) ' , snd8 (2:8) 
writadl, * ) f *t8 , 'D(FLUID) ' . snd8(2:8) 
alsaif (fst8.aq. ' VECTO'lthen 

writadl,*)' DISPLACEMENT (SORT1, PUNCH) = ALL’ 

else 

writadl. 1121) fst8. sndS. trd8, frtB, f the , aix8 , sth8 , 

4 athS . nth8 , tch8 

and if 
goto 100 
200 closedO) 
close (11 ) 

1121 format (10a8) 
stop 
and 

Program pch2res.f 

C1234412 34 5 6789012 

program pch2ra# 

character *72 title, subtitle. label, analysis type, datatype 

character* 30 ifile.ofile 

character*15 subcase, eigen 

character*8 lease, inumbar , sn, mn 

character *6 mode, con t 

write (*,*) 'ENTER PUNCH FILE' 

raad(*, l)ifile 

write (*,*> 'ENTER NUMBER OF DATA POINTS' 
read ( * , * ) numdata 

open (uni t=ll, f ile* if lie, status* ' old' ) 

1101 read(ll,2,and=205) title, iline 
read (11, 2) subtitle, iline 
read (11, 2) label. iline 
read (11, 2) analysis type, iline 
read (11, 2) datatype, iline 

1 format (a) 

2 format (a72 , 18) 
c 

raad (11, 3) subcase, icase, iline 

raad (11, 4 1 eigen, freq.mode, inumbar, iline 

3 format (al5, 5x, a8, 44x, i8) 

4 f ormat (al5, E13 . 4 , 2x, a6 , a8 , 28x, i 8 ) 

if (analysistype (2:11) .eq. ’ EIGENVECTO ' ) then 
c extract job name 
i=0 

105 i=i + l 

if(ifile(i:i) .ne. ’ . ' . and. ilen . It . 20) then 


of ile (1: i ) sif ile (1 : i) 

ilen=i 

goto 105 

else 

endif 

c extract mode number 
i=0 
im=0 

106 i=i+l 

i f I inumbar ( i : i ) . ne . ' ' . and . im . le . 8 ) then 
im=im+l 

mn (im: im) = inumber (i : i ) 
goto 106 

alsaif (i. It. 8) than 
goto 106 

else 

endif 

c extract subcase number 
i=0 
Ib=0 

107 1=1+1 

if (lease (id) .ne. ' ' .and.is.le. 8) then 
is=is+l 

sn (is: is) * lease ( i : i) 
goto 107 

alsaif (i . It . 8) then 
goto 107 

else 

endif 

c create output file name 
write ( * , * ) inumbar , mn 
write ( * , * ) lease , sn 

of ila= i filed : ilan) // '_roode' //mn(l:im) / / ’ .dis. ' //sn(l:is) 

else 

write ( *, *) analysistype (1:11) 
of ile= ' error . dis . 1 ' 
endif 

c 

defmax=l 

nwidth*6 

ndmax= int ( numdata /2 ) 

opan (unit= 12 , file»o£ila, status* 'unknown' ) 
write (12, 1003) eigen, freq 

write(12, 1 111) numdata, numdata, da f max, ndmax, nwidth 
write (12. 5) title 
write (12 , 5) subtitle 
5 format(a72) 

1003 format (al5,E13. 4) 
do 100 1=1 , numdata 

read (11 , 1001 ) cont , node, typevar , al, a2 , a3 , iline 
readdl, 1002) cont, a4,a5,a6, iline 
write (12, 1112 Inode, al,a2, a3, a4, a5, a6 

1111 f ormat (219, el5. 9, 219) 

1112 format (18, (5al3.7) ) 

1001 format (a6, 18, a4, 3 (5x,el3.4) ,18) 

1002 format (a6, 12x, 3 (5x,el3.4) . i8) 

100 continue 

close (12) 
goto 1101 
205 close(ll) 
stop 
end 


Program convert. f 

program convert 
c 

c Program to convert xxx.pch file from NASTRAN forced- response 
c analysis to PATRAN- readable XY-Plot files or TECPLOT 
c formatted data files, 
c 

c Input consists of punch file from NASTRAN forced-response 
c analysis. The total number of nodes (maxn) in the model, the 
c number of nodes calculated (iwant), the number of frequencies 
c (roaxfreq) and the output format must be programed prior to 
c compiling, 
c 

c Output consists of XY-Plot or TECPLOT date for displacement 
c magnitudes at a point vs. frequency. XYPlot files contain 
c magnitude or phase information for all nodes. Tecplot files 
c contain infromation for phase and magnitude for only one node 
c 

c Written by C.M.Femholz 48221 
c 

c Declarations 
c 

implicit real*8(a-h,o-y> 
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implicit complex(z) 
c 

parameter (maxn=5643 , maxf req=126 , iwant=3 , pi= 3 . 141592654 J 
c 

character*72 title, subtitle, label, analysis type, datatype 
character* 30 if ile, of ileRe, ofilelm, ofileu,ofilev, ofilew, 
k ofilex,of iley, ofilez 

character*15 subcase, point , displace 
character* 8 i case , inumber , sn, mn 
character* 6 mode,cont 
character*4 typevar 
c 

dimension disr (6) , disi (6) , f requency (maxfreq) , node(maxn) , 
& freqRe {maxf req, maxn , 6) , f reqlm {maxf req, maxn, 6) , 

St udata {maxn, maxf req) , utheta {maxn, maxf req) , 

k vdat a (maxn, maxfreq) , vtheta (maxn , maxf req) , 

St wdata {maxn, maxfreq) , wtheta (maxn, maxfreq) , 

St uphase {maxn, maxf req) , vphase (maxn, maxfreq) , 

St wphase {maxn, maxf req) , size( 3) 

c 

logical patran, tecplot 
c 

c Begin program 

c 

c Choose output type 
c 

patran= . FALSE . 
teCplOtr.TRUE. 
c 

write{ *,*) ’ Enter punch file name' 
read {* , 1 ) if ile 
1 format(a) 


Open punch file 

open (uni t= 10, file=ifile,status=' old' ) 
do 50 i=l,iwant 


c Read headers from punch file 
c 

read (10, 1000) title, i line 
read{10, 1000) subtitle, iline 
read{10, 1000) label , iline 
read{10, 1000) analysistype , iline 
read (10, 1000 ) datatype , iline 
read(10, 1005) subcase, lease, iline 
read (10, 10 10) point, node (i) 
c 

c Read in frequencies, displacements 
c 

do 100 j=l , maxfreq 
c 

read (10, 1100) frequency (j) , typevar, 
k (disr(k) , k=l,3) , iline 

read (10, 1105)cont, (disr Ik) , k=4, 6) , iline 
read (10, 1105) cont, (disi (k) , k=l, 3) , iline 
read (10, 1105)cont, (disi Ik) , k=4, 6) , iline 
c 

do 9 k=l, 3 

zdisp=cmplx (disr (k) , disi (k) } 
partone=real (zdisp) 
parttwo=imag { zdisp ) 

size (k) =SQRT (partone*partone+parttwo ‘part two) 

9 continue 

c 

c Store real and imaginary displacements for each frequency in an 
c array. Note: used if PATRAN fringe plots of displacements for 

c a given frequency are desired, 
c 

do 160 k=l, 6 

freqRe ( j , i , k) =disr (k) 
f reqlm ( j , i , k) =disi (k) 

160 continue 

c 

c Store displacement magnitudes for each node 
c 

if (patran) then 

udata (i, j ) =size { 1) 
vdat a (i , j ) =size { 2 ) 
wdata (i, j ) =size ( 3 ) 
endif 
c 

if (tecplot) then 

udata (node ( i } , j ) =size (1 ) 
vdata (node (i) , j ) =size (2 ) 
wdata {node (i I , j ) =size (3 ) 
endif 


c 

c Store displacement phase angles for each node 
c 

if (disr (1) .EQ. 0.0) then 

if (disi (1) .GT. 0.0) then 


utheta (i, j) =90.0 
elseif (disi (1) .LT. 0.0) then 
utheta {i , j ) =-90 . 0 
elseif (disi (1) .EQ. 0.0) then 
utheta (i , j ) =0 . 0 
endif 
else 

utheta (i, j )= {180. 0/pi) *ATAN (disi (1) /disr (1) ) 
endif 

if (disr (2) . EQ .0.0} then 
if (disi<2) .GT.0.0) then 
vtheta { i, j } =90 . 0 
elseif (disi (2) .LT. 0.0) then 
vtheta ( i , j ) =- 9 0 . 0 
else 

vtheta (i , j ) =0 . 0 
endif 
else 

vtheta(i,j)= (180. 0/pi) *ATAN (disi (2) /disr (2) ) 
endif 


if (disr (3) .EQ. 0.0) then 
if (disi (3) .GT.0.0) then 
wtheta (i , j } =90 . 0 
elseif (disi (3) .LT.O .0} then 
wtheta (i, j)=-90.0 
else 

wtheta (i , j ) =0 . 0 
endif 
else 

wtheta ( i, j)= (180. 0/pi) *ATAN (disi (3) /disr (3) ) 
endif 
c 

if (tecplot) then 

uphase ( node (i) , j ) =u t he t a ( i , j ) 
vphase ( node ( i ) , j ) =vtheta ( i , j ) 
wphase (node(i ) , j ) = wtheta ( i , j ) 
endif 
c 

100 continue 

50 continue 
c 

close (10) 
c 

c Dummy print 

write (* , * ) title 
c 

c Magnitude output files for u,v, and w displacements 
c 

if (patran) then 

of ileu= ’magu.xyd' 
of ilev= ‘magv.xyd' 
o f i lew= ‘ magw . xyd ' 
elseif (tecplot) then 
of ileu=' magu.pl t ' 
o f i lev= ' magv .pit' 
of ilew=' magw. pit ' 
endif 
c 

c Phase angle ouput files for u,v, and w displacements 
c 

if (patran) then 

o f i lex= ' phau . xyd ' 
of iley= 'phav. xyd' 
of ilez= 'phaw.xyd' 
elseif (tecplot) then 
of ilex= 'phau. pit ' 
of iley= 'phav. pit ' 
of ilez= 'phaw.plt * 
endif 
c 

open (unit=40, file=ofileu, status= ' unknown ' ) 
open (unit=50, file=of ilev, status* ' unknown ' ) 
open <unit=60 , file=ofilew, status= ' unknown ' ) 
open (unit=70, f ile=of ilex, status= 'unknown’ ) 
open | unit =80, f ile=ofiley, status= ' unknown ' ) 
open (unit=90, t ile=ofilez, status= ' unknown' ) 
c 

if (patran) then 
do 300 1=1 , iwant 
c 

write{40, 1300) ' XYDATA , U-DISP, NODE' , node(l) 
write (50, 1300) ' XYDATA, V-DISP, NODE' ,node(l) 
write (60, 1300) 'XYDATA, W-D ISP, NODE ', node { 1 ) 
write (70,1300) ' XYDATA, U- PHASE NODE' , node (11 
write (80, 1300) 'XYDATA, V- PHASE NODE' ,node(l) 
write (90, 1300) 'XYDATA, W- PHASE NODE' ,node(l) 
c 

do 350 m=l, maxfreq 

write (40, 1305) frequency(m) ,udata(l,m) 
write (50, 1305) frequency(m) , vdata{l,m) 
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writ® (60, 13051 f raquancy (m) , wdata(l,m) 
writ® (70, 1305 ( frequency (m) ,utheta(l,m) 
writ® (80,1305) frequency (m) , vtheta ( 1 , m) 
write (90, 1305) frequency (m) , wtheta(l,m) 
350 continue 

300 continue 


c 


14 


do 14 n=40, 90, 10 

write(n, 1310) 'END' 
continue 


c 

c 


c 


11 

c 


12 

c 


2 


450 

c 

c 


13 


elseif (tecplot) then 


write (40, 1400} 'TITLE 
write (50, 1400} 'TITLE 
write (60, 1400) 'TITLE 
write (70, 1401) 'TITLE 
write (80, 1401) 'TITLE 
write { 90 , 1401 } ' TITLE 


•Magnitude for u- displacements ‘ ' 
•Magnitude for v-displacements* ' 
•Magnitude for w- displacements* ' 
•Phase Angle for u- displacements 
•Phase Angle for v-displacements 
•Phase Angle for w- displacements 


do 11 n=40 ,60,10 

write (n, 1405 )' VARIABLES = ■ Frequency , ‘Displacement ■ ' 
write (n, 1410 )' ZONE T= "Numeric Solution 1 , I =',maxfreq 
continue 


do 12 n=70, 90, 10 

write(n, 1406) 'VARIABLES = • Frequency* ,* Phase ■ ' 
write (n, 1411) 'ZONE T=‘ Phases*, I =’,maxfreq 
continue 


write {*,*} 'Enter desired node number for analysis' 
read (*,2)1 
format (i) 

do 450 m=l,maxfreq 

write (40, 1415) frequency (m) ,udata(l,m) 
wr i te ( 50 , 141 5 ) frequency (m) , vdata ( 1 , m ) 
write (60, 1415) frequency (m) , wdata(l,m) 
write(70, 1415) frequency(m) ,uphase (l,m> 
write (80, 1415) frequency(m) , vphase(l,m) 
write (90, 1415 ) frequency (m) , wphase ( 1 , m) 
continue 


endif 


do 13 n=40, 90, 10 
close (n) 
continue 


1000 format <a72,i8) 

1005 format (al5, 5x, a8, 44x, i8) 

1010 format (al3, 5x, i8, 46x,i8) 

1100 f ormat (4x, el3 . 4 , al , 3 ( 5x, el3 . 4) ,18) 
1105 format (a6, 12x, 3 (5x, el 3 .4) ,18) 

1200 format (al5,el3 .4) 

1205 format (2i9, el 5.9,219) 

1210 format (a72) 

1215 format (i8, (5el3 .7) /el3.7> 

1300 format (al9 , 14) 

1305 format (fi0.3,2x,fl0.6) 

1310 format (a3) 

1400 format(a39) 

1401 format (a41 ) 

1405 format(a38) 

1406 format (a31 ) 

1410 format (a30, i 5 ) 

1411 format (a20,i5) 

1415 format ( f 10 . 3, 2x, el 3 . 5) 

1420 format (a2 5, i5) 


c 

9999 stop 
end 


Program reader.f 


program reader 
c 

c Program to compare analytical and numerical displacement 
c calculations for cubic geometry, fluid only, 
c 

c Input files include: 

c NASTRAN- genera tad neutral file: used to determine 

c geometric location of each node in the model. Neutral 
c file should contain only node location information, 

c NASTRAN-generated punch file: used to determine the 

c displacement of each node, as calculated by NASTRAN. 
c 

c Output files: 

c read. out: includes information about the maximum 

c displacement error for each direction and the maximum 


c displacement error magnitude for the system, 
c error. dis.l: PATRAN-readable displacement file that 

c can be used to display the displacement errors for the 
c model in fringe plot form. 

c 

c It is necessary to specify the number of nodes (maxn) and the 
c number of elements (maxe) in the model prior to corq>iling. 
c 

c Declarations 
c 

implicit real*8(a-h,o-z) 
c 

parameter (maxn=4961, maxe=1000) 
c 

character*lS zfile.pchfile, eigen 
character* 8 count , typevar 
character za(maxn)*l 
c 

dimension x (maxn) ,y (maxn) , z (maxn) , 

4 ap rod (maxn) , 

4 px(maxn) ,py(maxn) ,pz (maxn) , 

4 node (maxn) , error (maxn) , 

4 ndf (maxn) , junk (6) 

c 

c Begin program 
c 

write<*, *) 'Enter neutral file name' 
read(*, 1) zfile 

write (*,*)' Enter punch file name' 
read { * , ljpchf ile 
1 format (a) 
c 

c Define geometry (cube) 
c 

side=5.0 
pi=3. 141592654 
idn=l .0 
c 

c Mode shape of interest <nnn=x-dir , mmm=y-dir , kkk=z-dir) 
c 

nnn=l 

kkk=-3 

c 

c Open neutral file 

open (uni t=10, f ile=zf ile, status= ' old' ) 
c Open punch file 

open (unit=20, f ile=pchf ile, status= ' old' ) 
c open output file 

open (unit=30, files ' read. out ' , status- ' unknown ' ) 
c 

c Sort loop to determine maximum displacement in punch file 
c Punch file should contain data only for mode of interest 
c 

bigx=0 

bigysO 

bigz=0 

c 

do 50 j=l , ( maxn- 1 ) 
c 

read (20, 1004 ) cont, node ( j ) , typevar, 

4 px( j) ,py(j) ,pz(j) ,iline 

read | 20, 1005) cont, a4, a5, a6, i line 

c 

if (ABS(px{j) ) .GT.ABS(bigx) ) bigx=px ( j ) 
if (ABS(py(j) ) .GT.ABS(bigy) ) bigy=py(j) 
if (ABS(pz(j) ) .GT.ABS(bigz) ) bigzspz ( j } 
c 

50 continue 
c 

if (bigx.EQ.0.0) bigx=1.0 
if (bigy.EQ.0.0) bigy=1.0 
if (bigz . EQ . 0 . 0) bigz=1.0 
c 

write (30, 1) 'Error: Analytic vs. Numeric Solutions' 

write(30, 1015) 'Cubic geometry, mode ' , nnn,mmm, kkk 
write (30, 1010) 'bigx = ',bigx 

write ( 30, 1010 )' bigy = ',bigy 

write(30, 1010} 'bigz = ' , bigz 

c 

rewind (unit =20) 
c 

100 read ( 10 , 1000 ) idpacket , idn , iv , kc , nl , n2 , n3 , n4 , n5 
c 

if (idpacket. eq. 99) then 
close (10) 
close 120} 
goto 500 
c 

elseif ( idpacket . eq. 1} then 
c 

c Determine location of grid point 
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read(10, 1001) x (idn) ,y( idn) , z (idn) 

read (10, 1002) icf , za(idn) , ndf ( idn) , ncnf ig, ncid, 

& ( junk(i) , i=l, 6) 

c 

c Determine analytic displacement of grid point, mode (1,1,1) 
ax=sin( (nnn*pi*x ( idn) ) /side) 
ay=sin( <mmm*pi*y ( idn) ) /side) 
az=sin( (kkk*pi*z (idn) ) /side) 
aprod( idn) =ax*ay*az 
c 

c Determne grid point displacement from punch file 
read(20, 1004} cont, node(idn) , typevar, 

& px ( idn) ,py (idn) , pz (idn) , iline 

read (20, 100S) cont, a4 , a5 , a6 , iline 
px (idn) =px (idn) /bigx 
py (idn) =py (idn) /bigy 
pz (idn) =pz (idn) /bigz 
c 

else 

write (*,*)' Invalid packet ID' 
stop 
endif 
c 

goto 100 
c 

c Dummy print 

500 write (*,*) 'write to read. out completed' 
c 

c Calculate displacement error at each node, max error for system 
c 

bigerror=0 

c 

do 250 m= 1 , maxn 
c 

error (m) =aprod (m) -px (m) 
c 

if ( ABS( error (m) ) .GT. bigerror } then 
bigerror=error (m) 
endif 
c 

250 continue 
c 

write ( 30 , 1010) ' Max error '.bigerror 
c 

close ( 30 ) 
c 

c Write PATRAN input file 
c 

eigen- ' $EIGENVALUE 
f req=bigerror 
defmax=l 
nwidth=6 

ndmax=int (maxn/ 2 ) 
c 

open (uni t=40, f ile= ' error . dis . 1 ' , status=' unknown' ) 
write (40, 1100) eigen, freq 

write (40, 1101) maxn, maxn, defmax, ndmax, nwidth 
write (40, 1102) ' STITLE GOES HERE' 
write (40, 1102) ' $SUBTITLE=LOAD_CASE_ONE ' 
c 

do 200 1=1, maxn 

write (40, 1103) node (1) , error (1) ,0.0, 0.0, 0.0, 0.0, 0.0 
200 continue 
c 

close (40) 
c 

1000 format (i2, 8i8) 

1001 format (3el6. 9) 

1002 format (il, lal, 3i8, 2x, 6il) 

1003 format (4el6.9) 

1004 format (a6, i8 , a4, 3 ( 5x, el 3 . 4 } , i8 ) 

1005 format (a6, 12x, 3 ( 5x, el3 . 4) , i8) 

1010 format (al0,el8. 9) 

1015 format (a21,i2, ix, i2, lx, i2) 

1100 format (al5,el3. 4) 

1101 format (2i9,el5. 9, 2i9) 

1102 format (a72) 

1103 format (iB, <5el3.7) ) 
c 

stop 

end 

Program set maker. f 

program setmaker 
c 

c Program to produce NASTRAN AC MOD L set cards for a cubic geometry 
c as well as fluid/solid grid set. If grid point is determined to 
c be a fluid point, grid co-ordinate remains "-l * . If grid point 
c is a solid, co-ordinate ID is changed to , 0.‘ (solid), 
c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


1 

2 

c 


c 

c 


c 

c 

c 

c 


20 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


Input consists of original GRID data only from the NASTRAN 
bulk data file. It is necessary to specifiy the total number 
of grids in the model (maxn) prior to compiling. 

Output consists of two files: 

grids. out: contains fluid/ solid grid points for .bdf file 

set555: contains ACMODL set cards. 

In the ACMODL output file, SET1=555 contains the solid grid 
points and SET1=666 contains the fluid grid points. Line 
continuation markers start at ■AAAAAAA". 

Assumptions: Structure nodes are assumed to lie in two planes 
only; one at z=0 and the other at z=5. If a fluid node as a z 
coordinate of either 0 or 5, it is assumed to lie on the fluid/ 
structure interface and will be included in the ACMODL card 
set. It is assumed that every structure node iB on the fluid/ 
structure interface. The structure nodes are assumed to be 
sequential (ie xxx THRU xxx) . 

Program written by C.M.Femholz (48221) 

Declarations 

implicit real*8 (a-h, o-z) 

parameter (maxn=5643) 

character* 8 grid, set, thru, cont,con2 
character*15 gridf ile, yf lie, zfile 

logical test , check 

dimension inode (maxn) ,m (8) 

Begin program 

write (*,*)' Enter GRID file name' 
read ( * , 2 ) gridf ile 

write {*,*)' Enter minimum structure grid point ID' 
read ( * , 1 ) nodemin 

write (*,*)' Enter maximum structure grid point ID' 
read ( * , 1 ) nodemax 
format (i ) 
format (a) 

yf ile= ' grids. out ' 
zf ile= ’ set555' 

Open file containing grid information 

open (uni t=10, f ile=gridf ile, status= ' old' } 

Open file for ACMODL gridset 

open(unit=20, file=zfile, status= ' unknown ' ) 

Open file for grids output 

open(unit=30, file=yf ile , status= ' unknown ' ) 

Read grid file, determine fluid and solid grids, write new gridset 

do 20 i= 1 , maxn 

read ( 10 , 1000 ) grid , node , xx , yy , z z , f s 
if ( node . LT . nodemin . OR . node . GT . nodemax ) then 
fs=-l 

write (30, 1000) grid, node, xx, yy, zz, fs 
el seif ( node. GE. nodemin. OR. node. LE. nodemax } then 
£s=0 

write ( 30, 1000 ) grid, node , xx, yy , zz , fs 

else 

write (*,*) 'WARNING: Grid type indeterminate' 

endif 
continue 

close ( 10] 
rewind! uni t=30) 

Write structure SET card 

isset=555 
set= ' SET1 ' 
thru= ' THRU' 

write(20, 1005) set, isset, nodemin, thru, nodemax 
Read grid file, determine which fluid points lie on interface 
jj=0 

do 50 i=l , maxn 

read (30, 1000) grid, node, xx, yy, zz, fs 
if ( node. LT. nodemin. OR. node. GT. nodemax) then 
if (zz.EQ. (0.0) .OR.ZZ.EQ. (5.0) ) then 
inode { j j ) =node 
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endif 

c 

endif 

c 

50 continue 
c 

close {30} 
c 

c Write fluid SET card 
c 

55 if set=666 

c 

c Logic for line continuation markers 
c Note: *A* = char(65), *+* * char(43) 

c 

m(l)=43 
do 60 i-2 , S 
m(i) =65 
60 continue 
c 

cont=char (m{l> } //char (m(2) ) //char |m{3}} //char (m (4) ) // 
a char (m(5) ) //char (m( 6) ) / /char (m (7) } / /char (m (8) ) 
con2=cont 
c 

c Write first seven fluid points to SET card (first line) 
write {20, 1010} set , if set , { inode { k } , k=l, 7 ) , cont 
c 

c Write groups of remaining fluid grid nodes eight at a time 
c 

imax=INT { ( jj-7) / 8 } 
i things 8 
icount=l 
checks .TRUE, 
c 

100 if (check) then 
contscon2 
tests. TRUE. 

SO if (test) then 

i=8 

if (m(i).GE.90) then 
tn(i)=65 

m(i-l)=ro(i-l) +1 
if (m(i-l) . LT. 90} then 
tests. FALSE. 

else 

isi-1 

if ( i . LT . 2 } then 

write (*,*)' Too many grid points’ 
goto 999 
endif 
endif 
else 

m(i) =m (i) +1 
tests. FALSE, 
endif 

85 goto 80 

endif 

c 

con2=char (m(l) ) //char <m(2) ) //char <m(3) ) //char fm|4) }// 
6 char (m{ 5) ) //char (ra( 6) ) //char (m (7) } //char {m(8} } 

c 

write (20, 1015) cont , ( inode (k) , k=i thing, ithing+7 J , con2 
c 

i things! thing +8 

icount=icount+l 

if ( icount .GT. imax) then 
checks . false . 
endif 
goto 100 
endif 
c 

c Write last line of file 
c 

cont=con2 

write (20, 1020) cont, (inode (k) , k=i thing, jj ) 
c 

close (20) 
c 

1000 format (a8, i8, 8x, 3F8.4, i 8 ) 

1005 format (a8, 2i8, a8, i8) 

1010 format (a8, 18, 7i6, a8) 

1015 format (a8, 8i8 , a8 ) 

1020 format (a8 , 8i 8) 
c 

999 stop 
end 

Program plate. f 

program plate 
c 


c Program to compare analytic and numeric displacement 
c calculations for the cubic fluid/ structure geometry, 
c Errors for w- displacement and pressure are determined, 
c 

c Input files include: 

c PATRAN-generated punch ( .pch} file. Numeric dis- 
c placements for each node are read from this file, 
c PATRAN-generated neutral (.out) file containing 
c node information only. Geometric node locations 
c are read from this file, 
c 

c Output files include: 

c error. out: contains information about the maximum 

c error occuring anywhere in the model . 
c error. dis.l: PATRAN- readable displacement file that 

c can be used to display the error results for the 
c model in fringe plot form, 
c 

c It is necessary to specify the number of nodes (maxn), 
c the number of elements (maxe) , and the mode shape of 
c interest prior to compiling, 
c 

c Assumptions: 

c Analytic model uses an uncoupled solution. NASTRAN 
c uses a coupled one. For the first structural mode, 
c coupling will occur between the two plates of the 
c model. The mode shape of the plate that is moving 
c will be superimposed (with some attenuation of the 
c amplitude) upon the other plate. To account for this, 
c displacement amplitudes in the two plates are handled 
c separately in this code. The rear plate (at z=0( is 
c assumed to be stationary. The error.dis file can be 
c modified to delete the nodes corresponding to the 
c rear plate, 
c 

c Written by C.H.Femholz (48221} 
c 

c Declarations 
c 

implicit real (a-h,o-z) 
c 

parameter (maxn=1573,maxe=1200,pl=3 .141592654) 
c 

character* 72 title, subtitle, label, analysistype, 

6 datatype, header 

character*30 pchfile, neuf ile, header2 
character* IS subcase, eigen 
character*8 lease, an, mn 
character*6 mode, cont 
character *5 solidtype, fluidtype 
character*! flustr 
c 

dimension px(maxn) , py ( maxn ) ,pz (maxn) ,node(maxn) , 
k w(maxn) , p (maxn) , serror (maxn) , ferror (maxn) , 

a nodofluid (maxn) , nodesolid (maxn) , 

St x (maxn) ,y (maxn) , z {maxn} , 

St junk(6| , itrash(9> 

c 

logical fluid 
c 

c Begin program 

c 

c choose mode shape of interest (m=x-dir , n=y-dir , k=z-dir) 
c 

m= -2 
n=2 
k=0 
c 

c Define cubic geometry 
c 

side=5. 0 
thick=0 . 0625 
c 

c FEM geometry 

c 

ii=0 

jj=0 

f luidtype= ' HEX8 ' 
solldtype= 'QUAD4 ' 
c 

write (*,*) 'Fluid or structure mode? (f/s) ' 
read (*,1) flustr 

write (*,*) 'Enter desired mode (eigenvalue} number' 
read ( * , 2 ) modenumber 

write!*, *) 'Enter punch {.pch} file name' 
read (*,1) pchfile 

write!*, *} 'Enter neutral (.out) file name' 
read (*,1} neuf lie 
1 format {a} 

c 

if ( flustr . EQ. ' f ' ) then 
f luid= .TRUE. 
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elseif { flustr . EQ. ' s ' } then 
fluid=. FALSE, 
endif 


c 

c Open punch file 

open (unit=10, f ile=pchfile, status= ' old' ) 
c Open neutral file 

open |unit=20, f ile=neufile, status= ’ old’ ) 
c Open error output data file 

open (units 30 , files' error . out ' , status= 'unknown ' > 
c 

write (*,*)' Enter min structure grid ID for front plate' 
read ( * , 2 l nodeminf 

write (*,*>' Enter max structure grid ID for front plate' 
read { * , 2 ) nodemaxf 

write {*,*} 'Enter min structure grid ID for back plate' 
read { * , 2 J nodeminb 

write {*,*}' Enter max structure grid ID for back plate’ 
read ( * , 2 } nodemaxb 
2 format (i) 


c 

c Read displacements from punch file 
c 


100 bigx=0 . 0 
bigy=0 . 0 
bigzf=0. 0 
bigzb=0.0 


read (10, 1000) title, iline 

read (10, 1000) subtitle, iline 

read{ 10, 1000) label , iline 

readtlO, 1000) analysistype, iline 

read(10, 1000) datatype, iline 

read|10, 1005) subcase, icase, iline 

read(10, 1010) eigen, freq.mode, inumber, iline 


c 

c 


do 150 j = l , maxn 

read ( 10, 1015) cont , node ( j) , typevar , px ( j ) , py ( j ) , pz ( j ) , 
& iline 

read (10, 1020) cont, a4, a5,a6, iline 


if (node( j ) .LT. nodeminb. OR. node < j) . LT . nodeminf . OR. 

& node ( j ) . GT . nodemaxb . OR . node < j ) . GT . nodemaxf ) then 

if (ABS(px| j) ) .GT.ABS(bigx) ) bigx=px ( j ) 
endif 


c 

if (ABS (py(j) ) .GT.ABS (bigy) ) bigy=py(j) 
c 

i f ( node ( j ) . LE . nodemaxb . AND . node { j ) . GE . nodeminb ) then 
if (ABS (pz ( j) ) .GT. ABS (bigzb) ) bigzb=pz(j) 
elseif (node ( j ). LE. nodemaxf . AND. node { j ) ,GE . nodeminf ) then 
if (ABS (pz ( j) ) .GT. ABS (bigzf ) ) bigzf=pz(j) 
endif 
c 

150 continue 


if ( inumber .NE. modenumber ) goto 100 
c 

close (10) 
c 

if (bigx. EQ. 0.0) bigx=1.0 
if (bigy . EQ . 0 . 0) bigy=1.0 
if (bigzf . EQ. 0. 0) bigzf=1.0 
if (bigzb. EQ. 0. 0) bigzb=1.0 
c 

c Write to error. out file 
c 

write{30, 1100) 'Analytic vs. Numeric Solutions' 
write{30, 1100) 'Cubic Fluid/Structure Geometry’ 
write{30, 1105)maxe, solidtype, ' , ' , fluidtype, 'elements' 
write{30, 1110)maxn, 'nodes' 


if (fluid) then 

write (30 , 1130) ' Fluid analysis' 
elseif (.NOT. fluid) then 

write ( 30, 1130) ' Solid analysis' 
endif 
c 

write { 30 , 111S )' bigx = ‘ , bigx 
write I 30, 1115) 'bigy = ' , bigy 
write (30, 1115) 'bizgf= '.bigzf 
write (30. 1115) 'bigzb= '.bigzb 
c 

c Normalize displacements 
c 

do 200 j = 1 , maxn 

px(j)=px(j) /ABS (bigx] 
py< j)=py ( j) /ABS (bigy) 

if (node! j) . LE. nodemaxb . AND. node (j ) .GE. nodeminb) then 
pz ( j ) =pz ( j ) /bigzb 

elseif (node ( j ) . LE. nodemaxf .AND. node ( j ) .GE . nodeminf) then 


pz ( j ) =pz ( j ) /bigzf 
endif 
200 continue 
c 

c Read grid point locations from neutral file 
c 

read (20, 1200) (itrash(i) ,i=l,9) 
read(20, 1205)header 
read (20, 1200 ) (itrash(i) ,i=l,9) 
read (20, 1210 ) header2 
c 

300 read (20, 1200 ) idpacket , idn , iv, kc , nl , n2 , n3 , n4 , n5 
c 

if ( idpacket . EQ. 99 ) then 
close (20 ) 
goto 350 
c 

elseif (idpacket . EQ. 1) then 
c 

read (20 , 1215 )x ( idn) , y (idn) , z ( idn) 

read (20, 1220)icf,za, ndf , ncnf ig, ncid, { junk ( i ) , i=l , 6) 
c 

if (idn. LE. nodemaxb. AND. idn. GE. nodeminb) then 
c 

ii=ii+l 

nodesolid (ii ) = idn 
w(ii) =0.0 
c 

serror (ii) =w(ii) -pz (idn) 
c 

elseif { idn . LE . nodemaxf . AND . idn. GE . nodeminf ) then 
c 

ii=ii+l 

nodesolid (ii ) =idn 
ax=sin{nr*pi*x(idn) /side) 
ay=sin(n*pi*y(idn) /side) 
w{ ii ) =ax*ay 
c 

serror (ii)=w(ii)-pz(idn) 
c 

else 


nodef luid{ j j ) =idn 
ax=sin (m*pi*x ( idn) /side) 
ay=sin (n*pi*y ( idn) /side) 
az=cos (k*pi*z ( idn) /side) 

P l j j ) =ax*ay*az 

f error ( j j ) =p{ j j ) -px(idn) 

endif 

else 

write ( * , * I 'WARNING: Invalid packet ID' 
endif 


goto 300 
c 

c Determine maximum errors in model 
c 

350 bigerrorb=0. 0 
bigerrorf=0. 0 
c 

do 400 i=l, ii 

if (nodesol id (i ) .LE. nodemaxb. AND. 

L nodesolid { i) .GE. nodeminb) then 

if (ABS { serror (i) ) .GT. ABS (bigerrorb) ) then 
bigerr orb= serror ( i } 
endif 

elseif (nodesolld(i) .LE. nodemaxf .AND. 

& nodesolid ( 1) .GE. nodeminf ) then 

if (ABS ( serror ( i) ) .GT. ABS (bigerrorf) ) then 
bigerrorf =serror ( i } 
endif 
endif 
400 continue 
c 

write (30, 1120} 'Maximum error in front plate ', bigerrorf 
write{30, 1120) 'Maximum error in rear plate ', bigerrorb 
c 

bigerror=0.0 

c 

do 450 i=l , j j 

if (ABS (f error (i) ) .GT. ABS (bigerror) ) then 
bigerror= f error ( i ) 
endif 
450 continue 
c 

write ( 30, 1125) ' Maximum error in fluid ', bigerror 
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close (30) 
c 

c open error. dis file 

open (unit=40, f ile= ' error . die . 1 * , status= 'unknown' ) 

c 

«igen= ' $ EIGENVALUE s' 

defmax=l 

nwidth= 6 

ndmax= I NT ( roaxn / 2 ) 

c 

if (.NOT. fluid) then 

subtitles' 5SUBTITLB = STRUCTURE ERROR' 
elseif (fluid) then 

subtitles ' SSUBT1TLE = FLUID ERROR' 
endif 
c 

write (40 , 1300) eigen, f req 

write (40, 1305 )maxn,maxn,def max, ndmax , n width 
write (40, 1310) title 
write (40, 1310) subtitle 
c 

if (.NOT. fluid) then 
c 

do 500 i=l,ii 

write (40, 1315) nodesolid (i) ,0.0, 0.0, 

L serror (i) , 0 . 0, 0 . 0, 0 . 0 

500 continue 

c 

elseif (fluid) then 
c 

do 550 i=l,jj 

write (40, 1315) nodefluid (i) , terror (i) , 

it 0.0, 0.0, 0.0, 0.0, 0.0 

550 continue 

c 

endif 

c 

close (40) 
c 

1000 format (a72, i 8 ) 

1005 fonnat(al5, 5x, a8, 44x, i8) 

1010 format (al5,el3 .4, 2x, a6, 18, 28x, 18} 

1015 format (a6, i8, a4, 3 (5x, el3 . 4) ,18) 

1020 format (a6, 12x, 3 (5x, el3.4) ,18) 

1100 format (a30) 

1105 format (i4 , lx, a5, al , lx, aS, lx, a8) 

1110 format (14, lx, a5) 

1115 format (a7 , f ) 

1120 format (a28, lx, f) 

1125 format (a22 , lx, f ) 

1130 format (al4) 

1200 format (12, 818) 

1205 format (a72 ) 

1210 format (a30) 

1215 format (3el6. 9) 

1220 format ( 11 , lal , 318 , 2x, 611) 

1300 format (al5,el3. 4) 

1305 format (219 , el 5 . 9,219) 

1310 format(a72) 

1315 format (18, (5el3.7) ) 
c 

stop 

end 

Program forres.f 

program forres 
c 

c Program to determine the analytic solution for the forced 
c frequency response of a cubic fluid/structure geometry, 
c Model consists of fluid cube with sides of length A having 
c structrual plates located at z = -A/2, A/2. 

c 

c No input files required. Fluid and structure material 
c properties must be programmed prior to conqpiling. 
c 

c Output consists of five files: 
c struc.eig: Eigen values of the structure 

c sf req.plt : Disp. vs. freq. of the structure at a point 

c sphas.plt: Phase angle vs. freq. of the struct at a point 

c pf req. pit: Pressure vs. freq. of a point in the fluid 

c pphas.plt: Phase angle vs. freq. of the fluid at a point 

c Output files are written in TEC PLOT format, 
c 

c A non-coupled model is assumed. Structure obeys non-homogeneous 
c plate equation of motion. Fluid obeys homogeneous wave equation, 
c 

c Written by C.M.Fernholz (48221) 
c 

c Declarations 
c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


1 

2 

c 


c 

c 

c 

c 

c 

c 


c 


implicit real (a-h,o-y) 
implicit complex (z) 

parameter (nfreq=1000, step=5. 0, pi =3 .141592654) 

dimension z<pnn(30, 30) , zdisp (nf req) , 
it zqpress ( 30, 30) , zpress (nfreq) 

logical once, print 

Begin program 

once*. TRUE, 
prints. TRUE. 

Fluid and structure material properties 
Two pi 

twopi=2 . 0*pi 

Length of a side (inches) 

A=5 . 0 

Thickness of plates (inches) 
hs0.0625 

Young's Modulus for the plates (psi) 

E=10.3e6 

Density of the plate (slugs/in**3) 
rhos=2 . 5383e-4 
Damping coefficient 
eta=0. 005 

Amplitude of input force (lbs) 

F=5 . 0 

Poission ratio for plate 
uus0.334 

Speed of sound in fluid (in/sec) 
co=13620. 0 

Density of the fluid (slugs/in**3 } 
rhof =1 . 17e-7 

Structure stiffness 

d= (E*h**3) / (12* (l-uu*uu) ) 

Point of interest 

write (*, 2) 'Enter x co-ordinate' 
read ( * , 1 ) x 

write (*,*)' Enter y co-ordinate’ 
read ( * , 1 ) y 

write (*, 2 )' Enter z co-ordinate' 
read ( * , 1 ) ez 

write (*,*}' Entered: ',x,y,ez 
format ( f4. 2) 
format (a) 

freq=0 .0 
skip=step* twopi 
zo= (0.0, 0.0) 
zi=(0. 0,1.0) 

open (unit=10, file=' struct .eig' , status= ' unknown ' ) 
open (unit=20, f ile=' sf req.plt ' , status= * unknown' ) 
open (uni t=25, f ile=' sphas .pit ' , status= ' unknown ' ) 
open (unit= 30, f il®='pf req.plt' , status= ' unknown ’ ) 
open (unit=35, f ile='pphas .pit ' , status- 1 unknown ' ) 

do 100 i=l, nfreq 

freq=i*skip 
zdisp(i)=zo 
zpress (i) =zo 

do 200 m=l , 20 
do 250 n=l , 20 

wo= ( (pi*pi* (n*n+m*m) ) / ( A*A) ) *SQRT(d/ (rhos*h> ) 
if (print) then 

write (10, 1000} 'Mode' ,m, ’ , ’ ,n, 'and 
4 Frequency' , wo /twopi 

endif 

zdenom=cmplx(wo*wo-freq*freq, 2. 0*eta*wo*f req) 
zqmn <m,n)=4.0*F*sin (m*pi/2 . 0 ) *sin (n*pi/2 . 0) / 

& ( A*A*rhos*h*zdenom) 

zdisp ( i ) = zdisp (i) +zqmn (m, n) * 
fc sin (m*pi*x/A) *sin !n*pi*y/A) 

alf=- ( f req/co) * { freq/co) ♦ 

4 (n*pi/A) * (n*pi/A| ♦ 

& (m*pi/A) * (m*pi/A) 

if (alf.LT.0.0) then 
alf =SQRT (ABS (alf ) ) 
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zalpha=cmplx(0.0, alf) 
else 

alf=SQRT (alf) 
zalpha=cmplx (alf, 0-0) 
endif 


zqpreBs {m, n) =- ( rhof *f req* f req*zqmn (m, n) ) /zalpha 
zpress ( i) =zpress ( i > +zqpress (m, n) * 

& sin (m*pi*x/A) * 

i sin (n*pi*y/A) * 

St (csin ( zalpha*ez} + 

& ( 1+ccos (zalpha*A) ) *ccos (zalpha*ez) / 

St csin( zalpha*A) ) 

250 continue 

200 continue 


c 

call tecplot (freq, zpress ( i) , zdisp(i) , 

St frequency, pmag, ppha.dmag, dp ha) 

c 

if (once) then 

write{20, 2000) 'TITLE = ‘Analytic W-Displacements * ' 
write{25, 2000) 'TITLE = ‘Analytic W-Phase Angles “ 
write{30, 2000) 'TITLE = ‘Analytic U-Displacements ‘ * 
write{35, 2000) 'TITLE = ‘Analytic U-Phase Angles •' 
c 

write (20, 2 005) 'VARIABLES = • Frequency *,“ Magnitude * ' 
write (2 5, 2010) 'VARIABLES = * Frequency* , ‘ Phase ■ ' 
write (30, 2005) 'VARIABLES = ‘ Frequency* , ‘Magnitude ■ ' 
write (35, 2010) 'VARIABLES = * Frequency \* Phase • ' 
c 

do 300 j=2Q ,35,5 

write { j , 2015) ' ZONE T= ‘Analytic Solution*' 

300 continue 

c 

once= . FALSE, 
endif 
c 

write (20, 2020 ) frequency, dinag 
write (25, 2020) frequency , dpha 
write (30, 2020) frequency , pmag 
write (35, 2020) frequency , ppha 
c 

print= . FALSE . 

100 continue 
c 

close ( 10 ) 
do 305 j=20, 35, 5 
close ( j ) 

305 continue 


1000 format (a4 , lx, i2 , al , i2 , lx, al3 , lx, fl0.2) 
2000 format (a34) 

2005 format (a35) 

2010 format (a31) 

2015 format (a26) 

2020 format ( f 10 .2, lx, el3 . 5) 


stop 

end 


subroutine tecplot (freq, zpres, zdisp, fre, pmag, ppha, dmag, c^pha) 
c 

c Subroutine to format data for use by TECPLOT 
c 

c Declarations 
c 

implicit real (a-h,o-y) 
implicit complex (z) 
c 

parameter (pi=3 . 14 1 592654 ) 
c 

c Begin subroutine 
c 

c Frequency output 

c 

f re=f req/ (2*pi ) 
c 

c Pressure output 
c 

preal=real (zpres ) 
pimag=imag (zpres ) 

pmag= SQRT (preal *preal+pimag*pimag) 

c 

if (preal . EQ . 0 . 0 ) then 
if (pimag. GT. 0. 0 ) then 
ppha= (pi/2.0) 

elseif (pimag . LT. 0 . 0) then 
ppha=- (pi/2. 0) 

else 

ppha=0 . 0 


endif 

else 

ppha=ATAN (pimag /preal ) 
endif 
c 

c Displacement output 
c 

dreal=real (zdisp) 
dimag= imag ( zdisp ) 

dmag=SQRT (dreal *dreal+dimag*dimag ) 
c 

if (dreal . EQ. 0. 0) then 
if (dimag. GT. 0.0) then 
dpha= (pi/2 .0} 

elseif (dimag. LT. 0.0) then 
dpha=- (pi/2 .0) 
else 

dpha=Q . 0 
endif 
else 

dp ha = AT AN ( dimag /dr ea 1 ) 
endif 
c 

return 

end 

c 


complex function ccos(z) 
c 

c Function to return the cosine of a coir^lex argument 
c 

c Declarations 
c 

implicit real (a-h,o-y) 
implicit complex (z) 
c 

c Begin function 
c 

argr=real (z) 
argi=imag ( z) 
c 

partone=cos (argr) *cosh (argi ) 
parttwo=-l*sin(argr) *sinh(argi) 
c 

ccos=cmplx (partone, parttwo) 
c 

return 

end 


complex function csin(z) 
c 

c Function to return the sine of a conplex argument 
c 

c Declarations 
c 

real argr , argi , partone, parttwo 
complex z 
c 

c Begin function 
c 

argr=real (z) 
argi=imag(z) 
c 

partone=sin(argr) *cosh(argi > 
par ttwo= cos (argr) *sinh(argi) 
c 

csin=cmplx (partone, parttwo) 
c 

return 

end 

Program compare. f 

program compare 
c 

c Program to compare analytical and numerical displacement 
c calculations of first normal mode for a cylinder- shaped 
c geometry. 

c 

c Input files: 

c PATRAN- genera ted neutral file: used to determine geometric 

c location of each node in the model. Neutral file should 

c contain only node location information, 
c NASTRAN-generated punch file: used to determine the 

c displacement of each node, as calculated by NASTRAN . 

c It is necessary to specify the total number of nodes (maxn) 
c and the number of elements (maxe) in the model prior to 
c compiling, 
c Output files: 

c error. out: includes information about the maximum 
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c displacement 

c error for the system. 

c error. dis.l: PATRAN- readable displacement file that can be 

c used to display the displacement errors for the model in 

c fringe plot form, 

c 

c Written by C.M.Fernholz (48221) 
c 

c Declarations 

implicit real*B <a-h, o-z) 
c 

parameter (maxn=6609, maxe=2440, pi=3 . 141592654 ) 
c 

character *15 zf ile,pchfile, eigen 
character *8 count , typevar 
character za(maxn)*l 
c 

dimension x (maxn) ,y (maxn) , z (maxn) , 
k ax (maxn) , ay (maxn) , az (maxn) , 

k px(maxn) ,py (maxn) ,pz (maxn) , 

k node (maxn) , error (maxn) , aprod (maxn) , 

k ndf(maxn), junk (6) 

c 

c Begin Program 
c 

write(* ,*) 'Enter Neutral File Name' 
read ( * , 1 ) zfile 

write (*,*)’ Enter Punch File Name' 
read( *, l)pchf ile 
1 format (a) 
c 

c Define Geometry 
axis=5 . 0 
radius=l . 0 


c 

c 

c 

c 

c 


c 

c 

c 


Zero-order Bessel function, first root 
root=2. 40482556 

Open neutral file 

open (unit=10, file=zfile, status='old' ) 

Open punch file 

open(unit=20, f ile=pchf ile, statues ' old' ) 

Open error data output file 

open (uni t=30 , files * error. out ’ , statues’ unknown' ) 

Sort loop to determine maximum displacement in punch file 


bigx=0 . 0 
bigy=0 . 0 
bigz=0 . 0 

do 50 j= 1 , maxn 

read (20, 1004 ) cont , node ( j ) , typevar, 

* pxl j) ,py( j) ,pz ( j) , iline 

read (20, 1005) cont, a4, a5, a6, iline 

1004 format |a6, 18, a4, 3 (5x, el3. 4 ) ,18} 

1005 format (a6, 12x, 3 ( 5x, el 3 . 4) , 18) 

if (ABS{px{ j) ) .GT.ABS(bigx) ) bigx=px(j) 
if (ABS(py(j) ) .GT.ABS(bigy) ) bigy=py { j ) 
if (ABS (pz ( j ) ) . GT. ABS (bigz ) ) bigz=pz(j) 

50 continue 


if (bigx. EQ. 0.0) bigx=1.0 
if (bigy.EQ.0.0) bigy=1.0 
if (bigz. EQ. 0.0) bigz=l . 0 
c 

c Dummy print 

write (*,*) 'Displacements read’ 
c 

c Read in data from neutral file 
c 

100 read { 10 , 1000 ) idpacket , idn , iv , kc , nl , n2 , n3 , n4 , n5 

1000 format (12, 818) 
c 

if ( idpacket. EQ. 99) then 
close (10} 
close (20} 
goto 500 
c 

elseif (idpacket . EQ. 1) then 
c 

c Determine location of grid point 

read (10, 1001) x (idn) ,y (idn) , z(idn) 
read(10, 1002) icf , za(idn) ,ndf (idn) , ncnfig,ncid, 
k { junk(i) , i=l, 6) 

1001 format (3el6. 9) 

1002 format (il, lal, 3i8, 2x, 6il) 
c 

c Determine analytical displacement of grid point 


rad= SQRT ( x { idn ) *x(idn) +y(idn) *y(idn) ) 
axial=sin ( (pi»z (idn) ) /axis) 
radial=BESSJ0 ( (rad*root) /radius) 
aprod { idn I =axial*radial 
c 

c Normalize grid point displacements from punch file 
px(idnl=px(idn) /bigx 
py(idn) =py (idn) /bigy 
pz(idn)=pz(idn) /bigz 
c 

else 

write (*,*)' Invalid packet ID’ 
stop 
endif 
c 

goto 100 
c 

c Calculate displacement error at each node, max error for system 
c 

500 bigerrorsO.O 
c 

do 250 rn= 1 , maxn 

error (m) =aprod (m) -px (m) 

if (ABS (error (m) ) .GT. ABS (bigerror) ) bigerror=error (m) 

250 continue 
c 

write (30, 1 )' Error : Analytic vs. Numeric Solutions' 
write ( 30, 1 ) ’Cylinder , fluid only, mode 1,0,1' 
write (30, *) ’bigx = '.bigx 
write (30, * ) 'bigy = '.bigy 
write (30, *} ’bigz = ',bigz 
write < 30, *) 'bigerror= '.bigerror 
c 

c Dummy print 

write (*,*) 'Errors determined' 
c 

c Write PATRAN input file 

c 

eigen= ' $ EIGENVALUE = ’ 

freq=bigerror 

defmax=l 

nwidth=6 

ndmax= I NT ( maxn / 2 ) 
c 

open (uni t=40 , files ’ error. dis.l' , status= ' unknown ' } 
write(40 , 1100) eigen, freq 

write (40, 1101 ) maxn, maxn, def max, ndmax.nwidth 
write (40, 1102) ’ STITLE GOES HERE' 
write (40, 1102) ' $ SUBTITLE^ LOAD_C AS E_ONE ' 
c 

do 200 1=1 , maxn 

write (40, 1103 )node (1 ) , error (1) ,0.0, 0.0, 0.0, 0.0, 0.0 
200 continue 
c 

close (40) 
close ( 30) 
c 

1100 format (al5,el3 .4} 

1101 format (2i9,el5.9,2i9) 

1102 format(a72) 

1103 format (18 , (5el3.7) } 
c 

999 stop 
end 

Program card.f 

program card 
c 

c Program to products NASTRAN ACMODL set cards for a cylindrical 
c geometry as well as fluid/solid grid set. If grid point is 
c determined to be a fluid point, grid co-ordinate ID remains 
c *-l*. If grid point is a solid, co-ordinate ID is changed to 
c '0* (solid) . 
c 

c Input data consists of original GRID data only from the 
c NASTRAN bulk data file. 

c 

c Output files: 

c ‘grids. out* contains fluid/solid grid points 

c ‘set555* contains ACMODL set cards. SET1 555 = solid grids, 

c SET1 666 = fluid grids on fluid/solid interface. Line 

c continuation markers start with ’AAAAAAA* 

c 

c Assumptions: program matches each structure grid point with its 
c closest fluid grid point. Model can only represent a fluid 
c volume surrounded by a structural shell. The shell only touches 
c the fluid on one side. 2D elements must have nodal configura- 
c tions compatible with 3D elements (eg QUAD4 2D elements with HEX8 
c 3D elements ) . 

c Program verifies that no fluid grid appears twice in the ACMODL 


78 



c fluid set card. Structure grid points are assumed to be 
c sequential lie xxx THRU xxxj . 
c 

c Written by: Christian M. Fernholz 4-8221 

c 

c Declarations 

c Note: necessary to assign maxn=number of grids in model, 


implicit real *8 { a-h, o-z) 
c 

parameter (maxn=868} 
c 

character * 8 grid , set , thru , cont , con2 
character*15 gridf ile,yf ile, zf ile 
character*! type (maxn) 
c 

logical repeat , test , check 
c 

dimension node (maxn) ,m (8) , 
k yf {maxn} , xf (maxn) , zf (maxn) , nodf (maxn) , 

k ys (maxn) , xs (maxn) , zs (maxnl , nods (maxn) , 

Sc y tempi (maxn) , xtempl (maxn) , ntempl (maxn) , 

Sc xtemp2 (maxn) , ntemp2 (maxn) 

c 

c Begin program 
c 

write (*,*) 'Enter grid file name' 
read( * , 2) gridf ile 

write (*,*) 'Enter minimum structure grid point ID' 
read ( * , 1 ) nodemin 

write (*,*) 'Enter maximum structure grid point ID' 
read ( * , 1 ) nodemax 

1 format (i) 

2 format (a) 

c 

zf ile= ' set555 ' 
yf ile= 'grids . out ' 
c 

pi=3. 141592654 
rad= 1 . 0 
c 

c Error tolerances 
toly=0 . 001 
tolz=0 . 001 
tolr=0 . 001 
c 

c Initialize node array 
do 10 i = l , maxn 
node (i) =9999 
10 continue 


c Open file containing grid information 

open <unit=10, f ile=gridf ile, status= ' old' ) 
c Open file for AC MOD L set cards 

open (uni t=20 , file=zfile. status^ 'unknown' ) 
c Open file for solid/ fluid grid output 

open (unit=30, file=yf ile, status= 'unknown' ) 
c 

c Read grid file, determine fluid/solid grids, write new gridset 

c 

do 20 i=l,maxn 

c 

read (10, 1000} grid, nod, xx, yy, zz, fs 
c 

if (nod. LT. nodemin. OR. nod. GT. nodemax) then 
f s = -l 

write ( 30, 1000) grid, nod, xx, yy, zz, fs 
elseif ( nod. GE . nodemin. OR. nod. LE. nodemax! then 
fs = 0 

write (30, 1000 ) grid, nod, xx, yy, zz, fs 

else 

write( *,*) 'WARNING: Grid type indeterminate' 

end if 
20 continue 
c 

close { 10 ) 
rewind { uni t=30) 
c 

c Write structure grid point AC MOD L. set card 

c 

isset=55S 
set= ' SET1 ' 
thru= ' THRU' 

write (20, 1005) set, isset , nodemin, thru, nodemax 
c 

c Read grid co-ordinates (x,y,z) from new grid file 
c 

jrO 

k=0 

irad=0 

c 

do 50 i=l , maxn 


read(30, 1000) grid, nod, xx, yy, zz, fs 
i f ( nod . LT . nodemin . OR . nod . GT . nodemax ) then 

j=j+l 

xf (j)=xx 

yf < j)=yy 

zf ( j )=zz 
nodf ( j ) =nod 

radius= SQRT ( xx*xx+yy *yy ) 

if ( radius. LT. (rad+tolr) .AND. radius. GT. (rad- tolr) ) then 
irad=irad+l 
endif 

elseif ( nod. GE. nodemin. OR. n od, LE. nodemax J then 
k=k+l 
xs ( k ) = xx 
ys ( k) =yy 
zs (k) =zz 
nods ( k ) =nod 
else 

write (*,*) 'ERROR: Grid type indeterminate’ 

goto 999 
endif 
50 continue 


c Dummy print 

write (*,*) 'Grids successfully read' 
write {*,*)’ Fluid grids total: ',j 
write (*,*)' Structure grids total: ,k 

c Dummy print 

write (*,*) 'Number of surface structure nodes: ',k 

write (*,*)’ Number of surface fluid nodes: , irad 

c 

c Verify that two fluid grid points do not occupy same location 
c 

repeat= . FALSE. 
n=0 
c 

do 700 1=1, (j-1) 

do 710 ii= (i + 1 ) , j 
c 

if (Xf (i) .EQ.xf (ii) .AND. 

& yf (i) .EQ.yf lii) .AND. 

Sc zf (i) .EQ.zf (ii) ) then 

repeat= .TRUE . 
n=n+l 
endif 
c 

710 continue 
700 continue 
c 

if (repeat) then 

write (*,*) 'WARNING: ',n,' Fluid grids identical' 

endif 


c 

c Match structure grid points to closest fluid grid point 
c 

do 500 i=l , k 
c 

c Find all fluid grids with same z coordinate 
c 

11 = 1 

do 505 1=1, j 

if (Zf (1 > .GT. (zs (i) -tolz) . AND. 
k zf ( 1 > .LT. (zs { i ) +tolz> > then 

y tempi ( 11 ) =yf (1 ) 
xtempl ( 11 ) =xf (1 ) 
ntempl (ll)=nodf (1) 

11 = 11+1 

endif 

505 continue 
c 

c Find fluid grids with same y coordinate 
c 

do 510 mn=l , 11 

if (ytempl (mn> .LT. (ys(i)+toly) .AND. 

Ci ytempl (ran) .GT. (ys (i)-toly) ) then 

x temp 2 (mm)=xtempl (ran) 
n temp 2 (ram) = ntempl (ran) 
mm=mm+ 1 
endif 

510 continue 
c 

c Find fluid grid with closest x coordinate 

c 

err=1000000. 0 

do 515 n=l, mm 

er ror=ABS { xt emp2 ( n ) -xs ( i ) ) 
if (xtemp2 (n) .EQ.xs (i) ) then 
node ( i ) =ntemp2 ( n ) 
goto S00 


79 



elseif (error .LT. err) then 
node ( i ) =ntemp2 (n) 
err=error 
endif 

515 continue 
c 

500 continue 
c 

c Dummy print 

write(*, *} 'Structure /fluid grids matched' 
c 

c Check to see if any fluid points are repeated 
c 

repeat = . FALSE . 
n=Q 
c 

do 600 1=1, (K- 1 ) 

do 610 ii= < i + 1 ) , k 

i f (node ( i ) . EQ . node ( i i ) ) then 
repeat = .TRUE. 
n=n+ 1 
endif 

610 continue 
600 continue 


if (repeat) then 

write (*,*) 'WARNING: ', n, ' Nodes repeated in ACMODL card' 
endif 


c Write fluid grid point ACMODL set card 
c 

55 if set=666 

m ( 1 1 =43 
do 60 i=2 , 8 
m ( i ) =65 
60 continue 
c 

cont=char (m( 1 ) ) / /char (m (2) ) //char (m ( 3) ) //char (m (4 ) ) // 
fc char (m{5) >//char(m(6) ) //char (m(7) ) //char (m (8) ) 
con2=cont 
c 

c Write first line of fluid grid points 

write (20, 1010) set, if set, (node(n) , n*l, 7) ,cont 
c 

c Write remaining fluid grid points, eight at a time 

c 

imax=INT( Ik-7) /81 
c 

ithing=8 
i count =1 
checks .TRUE. 


100 if (check) then 
cont=con2 
tests .TRUE. 

80 if (test) then 

i = 8 

if (m(l).GE 90) then 
m ( i > >65 

mii-lj>m<i-l)*l 
if .LT.90} then 

tests FALSE 
else 

Ui-1 

if (i.LT.2) then 

write {•,*) 'Too many grid points' 
goto 999 
endi f 
endif 
else 

m ( i } =m ( i } ♦ 1 
tests . FALSE, 
endif 

85 goto 80 

endif 
c 

con2=char (m(l) ) //char (m(2) ) //char (m(3) ) //char (m<4) ) // 
fc char (m(5) ) //char (m(6) ) / /char (m (7) ) //char<m<8) ) 

c 

write(20, 1015)cont, (node(n) , n=ithing, ithing+7) ,con2 
c 

i t hi ng= i t hi ng 8 
icount=icount+l 
if I icount .GT. imax) then 
checks. FALSE, 
endif 
goto 100 
endif 


c Write last line of file 
cont=con2 

write (20, 1020)cont , (node(n) ,n=i thing, k) 


c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 


c 


c 


c 


c 

c 

c 

c 

c 


c 


c 


c 

c 

c 


c 


close! 30) 
close (20) 

1000 format (a8, i8, 8x, 3f 8 . 4, 18) 

1001 format (a8, i8. 8x, 3f8.4) 

1005 format (a8, 2i8 , a8, i8) 

1010 format (a8, 18, 718, a8) 

1015 format (a8, 818 , a8) 

1020 format (a8, 818) 

999 stop 
end 

Program modesm.f 

program modesnt 

Program to calculate natural frequencies of vibration for a 
cylindrical shell. 

No input files required. Shell material properties and 
geometric dimensions must be programmed prior to compiling. 

Output file freq. out contains mode numbers and natural 
frequency. 

Uses Epstein-Kennard theory and solution outlined in Leissa: 
Vibration of Shells, NASA SP-288 

Declarations 

implicit real (a-h,o-y) 
implicit complex (z) 
real 1 , nu, lam, kl, k2 , kO 

parameter (pi= 3 . 141592654) 

Begin program 

Cylindrical shell, dimensions, material properties (Al) 

Length (in) 

1=5.0 

Radius (in) 
a=1.0 

Thickness (in) 
h=0 . 0625 
Polssion's ratio 
nusO.334 

Young's Modulus (psi) 

E=10 . 3E6 

Density (slugs/in**3 ) 
rho=2 . 5383E-4 
Square root of -1 
zi=craplx(0 .0, 1. 0) 

open (uni t=10, files ' freq. out' , status= ’unknown' ) 
write (10, 1005} 'Cylindrical Shell Natural Frequencies' 
write ( 10, 1010 ) ‘ Epstein-Kennard Theory' 
writellO, 1015) 'Mode' , 'm' , 'n' , 'Roots (Re, im) : ' , 
t 'Root 1 (Hz)', 'Root 2 (Hz)', 'Root 3 (Hz) ' 


do 200 m=0,10 
do 300 n=0, 10 

if (m.EQ.O.AND.n.EQ.O) goto 300 

lam= (m*pi*a/l) 

Define Donnel-Mustari constants 

k2=l . + . 5* ( 3 . -nu) * <n*n+lam*lam) ♦ (h*h/ (12 . *a*a) } * 
fc (n*n+lam*lam> * (n*n+lam*lam) 

kl= .5* (1-nu) * ( (3 . +2. *nu) *lam*lam+n*n+ I n*n+lam* lam) * 
fc (n*n+lam*lam) + ( ( 3 .0-nu) / (1 . 0-nu) ) * (h*h/ (12. *a*a) ) * 

k (n*n+lam*lam) * (n*n+lam*lam) * (n*n+lam*lam) ) 

k0= .5* ( 1. -nu) * ( ( 1 . -nu*nu) * ( lam* lam* lam* lam) + 
fc (h*h/ (12. *a*a) ) * (n*n+lam* lain) * (n*n+lam* lam) * 

fc (n*n+lam*lam) * <n*n+lam*lam) ) 

Modifying constants for Epstein-Kennard theory 

delk2= (l+3*nu) / (1-nu) - (2-8*nu*nu+3*nu**3) ‘ lam*lam/ 
fc (2* (1-nu) * *2 ) - ( 19- 37 *nu + 19*nu*nu+nu**3 ) / 

fc (2* (1-nu) **2) -nu*nu* (n*n+lam*lam) / ( (1-nu) **2) 

delkl= (3+8*nu-5*nu*nu-nu**3) *lam*lam/ (2* (1-nu) ) + 
fc (2+nu) *n*n/2- (6+4*nu-8*nu*nu+3*nu**3)*lam**4/ 
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program £orced2 
c 

c Program to calculate natural frequencies of vibration for a 
c cylindrical shell and the forced response over a range of 
c frequencies, 
c 

c No input files required. NASTRAN-generated punch file is optional 
c if comparison to a numceric theory is desired. Shell material 
c properties and geometric dimensions must be programmed prior to 
c compiling, 
c 

c Output files: freq.out contains mode numbers and natural 

c frequencies. 

c displace. out contains displacement information 

c written in TEC PLOT form, 

c 

c Uses Donnell-Mustari theory and solution outlined in Leissa: 
c Vibration of Shells, NASA SP-288 
c 


Declarations 

implicit real (a-h,o-y> 
implicit complex (z) 
real 1 , nu, lam 

parameter (pi = 3 . 141592654 , maxf = 1000 , modes = 1 0 ) 
dimension u(maxf ) ,v(maxf } ,w(maxf ) ,p(maxf ) , freqfmaxf } 
logical print 
Begin program 

Cylindrical shell, dimensions, material properties (Al) 

Length (in) 

1=50.0 
Radius (in) 
a=10 . 0 

Thickness (in) 
h=Q . 0625 

Poission’s ratio 
nu=0 .334 

Young's Modulus (psi) 

E=10 . 3E6 

Density of the structure (slugs/in**3) 
rho=2 . 5383E-4 

Density of the fluid I slugs/in**3 ) 
rhof =1 . 170e-7 

Speed of sound in fluid (in /sec) 
co=13620 . 0 
Damping coefficient 
eta=0 . 005 

Applied force (lbs) 

Fo=- 50 . 0 

Square root of -1 
zi=cmplx (0. 0, 1.0) 

Common coefficient (wave speed in structure) 
cl2=E/ (a*a*rho* (l-nu*nu) } 

Frequency increment (Hz) 
step=5 . 0 

Starting frequency 

f r equency= 5 . 0 * (2 . 0*pi ) 

Print natural frequencies once 
print= . FALSE. 

Determine location of interest 

rr=a/2 
tt=0. 0 
yz=l/2 

format (a) 
format { f ) 

format (al 3 , f4 . 1, lx, f4. 1 . lx, f 4 . 1 ) 
tt=tt* (pi/180.0) 
if (print) then 

open! uni t=10, £ile=' freq.out' , status= unknown' ) 
write (10, 1005) 'Cylindrical Shell Natural Frequencies' 
write (10, 1010) ‘Donnell-Mustari Theory' 
write (10, 1015) 'Mode' , 'm' , 'n' , 'Roots (Re, Im) : ' , 
t 'Root 1 (Hz)', 'Root 2 (Hz)', 'Root 3 (Hz)' 

end if 

•Prep 1 tecplot output file 

open{unit=20, file= ' displace . out ’ , status= 'unknown' ) 
write(20,l) 'TITLE = ■ FORCED FLUID/ STRUCTURE CYLINDER RESP. - ' 
write(20,l) 'VARIABLES = * FREQ ’, 1 U * , ’ V ' , "W, * P ’ ’ 
write (20, 1) 'ZONE T= 'Analytic • ' 

do 100 k=l , maxf 

2U=cmplx {0 .0, 0 .0) 
zv=cmplx{0.0,0.0) 
zw=cmplx (0. 0, 0.0) 
zp=cmplx (0. 0, 0 .0) 
do 200 m=l, modes 
do 300 n=0, modes 

if (m.EQ.O . AND.n. EQ. 0) goto 300 

Calculate natural frequency for mode 

lam= (m*pi*a/l) 

Define cubic equation constants 

c2=l. + . 5* {3. -nu) * (n*n*lam*lam) ♦ (h*h/ (12 . *a*a) ) * 
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St (n*n+lam*lam) * (n*n+lam*lam) 

c 

els . 5* (1-nu) * { (3 . +2 . *nu) * lam*lam+n*n+ (n*n+lam*lam) * 

St (n*n+lam*lam) + { <3 . O-nu) / <1. O-nu) j * (h*h/ ( 12 . *a*a) ) * 

& (n*n+lam*lam) * (n*n+lam*lam) * (n*n+lam*lam) ) 

c 

cO= . 5* (1 . -nu) * ( (1 . -nu*nu) * (lam* lam* lam* lam) ♦ 

St (h*h/ (12 . *a*a) ) * (n*n+lam*lam) * (n*n+lam*lam) * 

k (n*n+lam*lam) * (n*n+lam* lam) ) 

c 

c Solve cubic aquation for omega aquarad. Usa mathod outlinad 
c in CRC Standard Mathematical Tables, 23rd ad. page 105 
c 

c= (1. 0/3.0) * (3.0 *cl-c2*c2) 

d= (1 . 0/27 . 0) * ( -2 . 0*c2 *c2*c2+9 . 0*c2*cl-27 . 0*c0) 
c 

deltas (d*d/4 .0) + (c*c*c/27.0) 
c 

if (delta. LT. 0.0} than 
par tones -d/2 . 0 
parttwosSQRT ( -1 .0*delta) 
zPP=cmplx (par tone, part two) 
zQQ=cmplx (par tone, -parttwo) 
else 

partone=-d/2. 0 
parttwo=SQRT (delta) 
zPP=cmplx{ (partone+parttwo) ,0.0) 
zQQ=cmplx( (partone-parttwo) ,0.0) 
andif 
c 

zPPszPP** (1. 0/3.0) 

ZQQ=Z QQ**(1. 0/3.0} 
c 

zf irst=- . 5* (zPP+zQQ) 
zsecnd=- . 5* (zPP-zQQ) *SQRT(3 .0) *zi 
c 

partonasraal (zf irst) + real (zsaend) 
parttwo=imag (zf irst) +imag (zsaend) 
parttre=real (zf irst) -real (zsaend) 
partfor=imag ( zf irst ) -imag (zsaend) 
c 

zrootl=zPP+zQQ+c2 / 3 . 0 
zroot2=cmplx ( (partona+c2/3 . 0) .parttwo) 
zroot3=cmplx ( (parttre+c2/3 . 0) , partfor) 
c 

zfreql=SQRT< zrootl *B/ (rho* (l-nu*nu) *a*a) ) 
zfreq2=SQRT(zroot2*E/ (rho* ll-nu*nu) *a*a) ) 
zfreq3=SQRT(zroot3*E/ ( rho* ( l-nu*nu) *a*a) ) 
c 

if (print) than 

write (10, 1000) 'Mode' , m, ' , ' ,n, 'and frequencies ' , 

St zfreql/ (2. *pi) ,zfreq2/ (2. *pi) , zfreq3 / (2 . *pi ) 

if (m.EQ. modes) close(lO) 
andif 
c 

if (imag (zfreql ) .NE. 0 . 0 .OR. 
fc imag(zfraq2) .NE.O.O.OR. 

fc imag (zfraq3) .NE. 0. 0) then 

write (*,*)' Fatal error: natural frequency is ' , 

A ‘imaginary’ 

goto 9999 
andif 


c 

call magnitude (zfreql , fraql ) 
call magnitude (zfreq2, freq2 ) 
call magnitude ( zfreq3 , freq3 ) 
c 

c Calculate force coefficient 


laro= (m*pi/l) 
c 

Pmn=a* (2 . 0/ (pi*l) ) *Fo*sin (m*pi/2 . 0) * (1 . 0+cos (n*pi) } 
c 

f reqAs-SQRT ( - (C12/ (2 . 0*a*a) ) * (nu-1 .0)* (nu*lam*lam* 

St a*a-n*n) /nu) 

freqB=-SQRT(- (C12/ (2.0*a*a) ) * (nu-1 .0) * (nu*lam*lam* 
t a*a+2 . 0*lam* lam*a*a*n*n) ) 

freqCl=SQRT( (C12/a*a) * ( lam*lam*a*a+n*n) ) 
freqC2=SQRTl- (C12 / (2.0*a*a) ) * ( lam*lam*a*a+n*n) * (nu-1.0) ) 
c 

c Calculate cofactors for zAmn, zBmn, zemn 
c 

zdet=cmplx ( ( frequency* frequency- freq3*freq3 ) , 

Si 2.0*eta*frequency*freq3) * 

Si cmplxl ( frequency * frequency- freq2*freq2) , 

Si 2.0*eta*f requency*freq2) * 

& cmplx( ( frequency ‘frequency- fraql* fraql) , 

St 2.0*eta*frequency*freql) 

c 

zcof A=cmplx { { frequency* frequency- freqA*freqA) , 

St 2 . 0*eta* frequency* freqA) 

zcofB=cmplx( ( frequency* frequency- freqB*freqB) , 


Si 2 . 0*et a* frequency* freqB) 

ZcofCscmplx( ( frequency*f requency- f reqCl*f reqCl ) , 

Si 2 .0*eta* frequency* freqCl) * 

Si cmplx( (frequency* frequency- f reqC2 * freqC2 ) , 

St 2 .0*eta*f requency* freqC2) 

c 

zAmn=- (Fmn/ (rho*h) ) * zcof A/ (zdet) 
zBmn= (Fmn/ (rho*h) )* zcof B/ (zdet) 
zCmn=- (Fmn/ (rho*h) ) *zcofC/ (zdet) 
c 

c Sum displacements 

c 

zu=zu+zAmn*cos (m*pi*yz/l) *cos (n*tt) 
zv=zv+zBmn*sinlm*pi*yz/l) *sin(n*tt) 
zwszw+zcmn*sin (m*pi*yz/l) *cos (n*tt) 

c 

c Solve for pressure coefficient 
c 

a lphsqr= ( frequency* frequency) / (co*co)- (m*pi/l) * lm*pi/l) 
c 

if (alphsqr . LT. 0. 0) then 
alph=SQRT ( -1 . 0 * alphsqr ) 
call cof imag (a, alph, n, bottom) 
zDmn= (zCmn*rhof ‘frequency* frequency) /bottom 
call bessimag(rr, alph, n, press) 
c 

elseif (alphsqr .GT. 0.0) then 
a lph= SQRT ( alphsqr ) 
call realcof (a, alph, n, bottom) 
zDmn= (zCmn*rhof *f requency* frequency) /bottom 
call bessreal (rr, alph, n, press) 
c 

elseif (alphsqr . EQ. 0 . 0) then 
write (*,*) 'Alpha = 0.0' 
c 

endif 

c 

c Sum pressure 
c 

zp=zp+zDmn ‘press *s in |m*pi*yz/l) *cos (n*tt) 
c 

300 continue 
200 continue 
c 

c Dummy print 

write (*, 4) 'Complex data calculated for ',k 
4 format | a28, i4) 

c 

c Convert complex displacements, pressure to magnitudes 

c 

call magnitude (zu,u(k) > 
call magnitude (zv, v(k> } 
call magnitude(zw,w(k> ) 
call magnitude (zp,p (k) ) 
c 

c Dummy print 

write ( * , 1105} zp 
c 

c Write to TBCPLOT file 
c 

freq(k) = frequency/ (2 . 0*pi) 
write{20, 1100) freq(k) , u (k) , v (k) , w (k> ,p (k) 
c 

c Increment frequency 

f requency= frequency + step* (2.0*pi) 
c 

100 continue 
c 

close (20) 
c 

1000 format (a4, lx, 12 , al , i2, lx, al5, lx, 3 ( f 12.4 , lx, f 3 . 1 , lx) ) 

1005 format (a37) 

1010 format (a22) 

1015 format (a4, 2x, al , 2x, al, lx, al4, 3 ( 5x, all) } 

1100 format ( f7 . 2, 4 < lx, el6 . 9) > 

1105 format (2{el6.9,lx) ) 
c 

9999 stop 
end 

c 

c *«******ft**«***********«4** a * aa ***»*»«*« a «* a * aaa#aaa##aa# « a## , 

subroutine magnitude (z, x) 
c 

c Subroutine to return the magnitude x of a complex argument z 
c 

implicit real (a-h.o-y) 
implicit complex (z) 
c 

partone=real ( z } 
parttwo*imag(z) 
c 

x=SQRT (partone*partone+part two ‘parttwo) 
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return 

end 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


c 


subroutine realcof ( a, alph, n, out ) 


Subroutine to return the denominator of the pressure coefficient 
term. Used when the argument of the n-th order bessel function 
is a real. 

implicit real (a-h,o-y) 
implicit complex (z) 

Square root of -1 
zi=cmplx(0.0,1.0) 

x=alph*a 

if ln.EQ.0) then 
terml=BESSJO (x) 
term2=BESSJl (x) 
elseif (n.eq.l) then 
terml=BESSJl (x> 
term2=BESSJ(2, x) 
else 

tennl^BESSJ ( n , x } 
term2=BESSJ { (n+1) ,x) 
endif 

Construct denominator 


out= (n/a) *terml-alph*term2 

return 

end 


c 

c 

c 

c 

c 

c 

c 


c 


subroutine cofimag (a, alpha, n, out) 


Subroutine to return the denominator of the pressure coefficient 
term. Used when the argument of the n-th order bessel function is 
complex. 


implicit real (a-h,o-y) 
implicit complex (z) 


x=alpha*a 

if (n.EQ.O) then 
terml=BESSlO (x) 
term2=BESSIl (x) 
elseif (n.EQ.l) then 
terml=BESSlO(x) 
term2=BESSI (2 , x) 
else 

terml=BESSI (n, xj 


term2=BESSI { (n+1 ) ,x) 
endif 
c 

c Construct denominator 
c 

out=alpha*term2+ (n/a) *terml 
c 

return 

end 

c 

c a************************************************************ 
subroutine bessreal (r, alpha, n, out) 
c 

c Subroutine to return the n-th order bessel function of a real 
c argument . 
c 

implicit real (a-h,o-y) 
implicit complex (z) 
c 

zi=cmplx (0 . 0, 1.0) 
x=r*alpha 
c 

if (n.EQ.O) then 
out=BESSJ0 (x) 
elseif (n.EQ.l) then 
outrBESSJl (x) 

else 

out=BESSJ (n, x) 
endif 
c 

return 

end 


c 

c •**********•«*«*««******••******«***«*******•**•*•**********•»• 
subroutine bessimag |r, alpha, n, out ) 
c 

c Subroutine to return the n-th order modified bessel function of 
c a real argument 
c 

implicit real (a-h,o-y) 
implicit complex (z) 
c 

x=r*alpha 

c 

if (n.EQ.O) then 
OUt=BESSI0(x) 
elseif (n.EQ.l) then 
out=BESSIl (x) 

else 

out=BESSI (n, x) 
endif 
c 

return 

end 
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Appendix C: NASTRAN File Useage Statistics 


File Size 
(MB) 

Analysis Type 

Calculated 

Eigenvalues 

Requested 

Frequencies 

Nodes 

DBALL 

SCRATCH 

4.686 

0.459 

SOL 103 

10 

- 

1331 

19.759 

5.947 

10 

- 

4961 

6.660 

0.975 

14 

- 

1573 

27.394 

8.004 

14 

- 

5644 

5.890 

0.500 

10 

- 

1449 

24.478 

6.554 

10 

- 

6609 

14.025 

3.228 

47 

- 

1953 

61.538 

17.564 

51 ! 

- 

8097 

23.400 

6.554 

SOL 108 

- 

250 

5644 

24.396 

7.406 

SOL 111 

20 

250 

5644 

48.931 

16.687 

SOL 108 

- 

250 

1953 


Figure 31: MSC/NASTRAN file useage data for DBALL and SCRATCH files. 
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